Compression of Remotely Sensed Astronomical Image Using Wavelet-Based Compressed Sensing in Deep Space Exploration

: Compression of remotely sensed astronomical images is an essential part of deep space exploration. This study proposes a wavelet-based compressed sensing (CS) algorithm for astronomical image compression in a miniaturized independent optical sensor system, which introduces a new framework for CS in the wavelet domain. The algorithm starts with a traditional 2D discrete wavelet transform (DWT), which provides frequency information of an image. The wavelet coefﬁcients are rearranged in a new structured manner determined by the parent–child relationship between the sub-bands. We design scanning modes based on the direction information of high-frequency sub-bands, and propose an optimized measurement matrix with a double allocation of measurement rate. Through a single measurement matrix, higher measurement rates can be simultaneously allocated to sparse vectors containing more information and coefﬁcients with higher energy in sparse vectors. The double allocation strategy can achieve better image sampling. At the decoding side, orthogonal matching pursuit (OMP) and inverse discrete wavelet transform (IDWT) are used to reconstruct the image. Experimental results on simulated image and remotely sensed astronomical images show that our algorithm can achieve high-quality reconstruction with a low measurement rate.


Introduction
Light-small is the development direction of the optical sensor system in deep space exploration [1][2][3]. This requires the system to execute multiple functions on the same hardware platform [2]. In our previous work, we proposed a miniaturized independent optical sensor system that can provide altitude and navigation information for spacecraft based on optical imaging measurement [3]. Those remotely sensed astronomical images are crucially important for future studies of celestial bodies [4]. Therefore, the system has the functions of image compression and transmission.
Remotely sensed astronomical images impose an urgent demand for compression considering the rapid increase of data, the complexity of the deep space environment, and the performance limitations of airborne equipment [5][6][7][8][9]. Remotely sensed images of extraterrestrial celestial bodies, especially planets, have some unusual features. A typical remotely sensed astronomical image is composed of background and a celestial body [10]. We can significantly reduce the amount of data in the background without losing important information such as geometric structure, edges, and surface textures, which are of great significance for navigation accuracy analysis and scientific exploration [11][12][13]. Therefore, we hope that an optical sensor system will compress remotely sensed astronomical images with a low compression ratio (CR) and low computational complexity while preserving as much information as possible.

Discrete Wavelet Transform
The wavelet transform is a powerful time-frequency analysis method that was developed to overcome the Fourier transform's lack of local capability in the time domain. It has the important property of good localization characteristics in time and frequency domains, and it can provide frequency information of each sub-band of a signal. The discrete wavelet transform (DWT) discretizes the scale and translation of the basic wavelet [31,32]. Figure 1 shows the realization of DWT and inverse discrete wavelet transform (IDWT). The image I 0 j−1 (or low-frequency sub-band LL j−1 of the j − 1 level) uses low-pass filter lk and high-pass filter hk to convolve each raw of the image, and performs downsampling ( 2 ↓ indicates column sampling, keeping all even columns). The columns of the image are convolved with lk and hk and downsampled ( ↓ 2 indicates row sampling, keeping all even rows). Finally, we obtain four sub-images with a size one-fourth of I 0 j−1 , including: Lowfrequency sub-band LL j , horizontal high-frequency sub-band HL j , vertical high-frequency sub-band LH j , and diagonal high-frequency sub-band HH j . The process of reconstruction is opposite that of decomposition. The sub-image is upsampled in the column direction ( ↑ 2 means row interpolation, i.e., 0 is filled between adjacent elements), convolved with low-pass filter Lk and high-pass filter Hk, and finally upsampled in the row direction and convolved with Lk and Hk [33]. This step is performed for each layer to obtain the reconstructed image I 0 j−1 .
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 15 elements), convolved with low-pass filter and high-pass filter , and finally upsampled in the row direction and convolved with and [33]. This step is performed for each layer to obtain the reconstructed image . The DWT of image ( , ) is [33] ( ) = √ ∑ ∑ ( , ) ( , ) , where is an arbitrary scale, × is the size of the image, ( ) is the approximate coefficient at scale , and ( , ) is the scaling function. DWT offers an effective method for 2D matrix transformation from the space domain to the time domain to separate signals according to low and high frequency. The low-frequency sub-band retains the content information of the original image, and the energy of the image is concentrated in the sub-band. The horizontal and vertical sub-bands maintain the high-frequency information in the horizontal and vertical directions, respectively. The diagonal sub-band maintains high-frequency information in the diagonal direction.

Compressed Sensing
CS shows [22,23] that if a signal is sparse in a certain transform domain, an observation matrix is used to project the signal linearly. By solving a convex optimization problem, the signal can be recovered [34]. CS effectively solves the limitation of the Nyquist theorem on the sampling rate and the waste of data caused by conventional compression.
This sparse representation of the signal is the basis of CS. Assuming that signal ∈ is sparse under the sparse basis = { , ⋯ }, the sparse representation is where is the length of . The matrix form is where is the sparse representation of . Τ is a measurement matrix of size x ( ≪ ). Sampling (measurement) is a process of linear projection. The compressed signal is where Θ is a sensing matrix. It is a problem of an ill-conditioned underdetermined equation. If Τ and satisfy the restricted isometry property (RIP), i.e., for any constant ∈ (0,1), then coefficients (the sparsity of is ( ≪ )) can be accurately reconstructed. The minimum value of is the RIP constant. The equivalent of RIP is that Τ and are orthogonal. The DWT of image I(x, y) is [33] WT ϕ (j) where j is an arbitrary scale, m × n is the size of the image, WT ϕ (j) is the approximate coefficient at scale j, and ϕ j (x, y) is the scaling function. DWT offers an effective method for 2D matrix transformation from the space domain to the time domain to separate signals according to low and high frequency. The low-frequency sub-band retains the content information of the original image, and the energy of the image is concentrated in the sub-band. The horizontal and vertical sub-bands maintain the high-frequency information in the horizontal and vertical directions, respectively. The diagonal sub-band maintains high-frequency information in the diagonal direction.

Compressed Sensing
CS shows [22,23] that if a signal is sparse in a certain transform domain, an observation matrix is used to project the signal linearly. By solving a convex optimization problem, the signal can be recovered [34]. CS effectively solves the limitation of the Nyquist theorem on the sampling rate and the waste of data caused by conventional compression.
This sparse representation of the signal is the basis of CS. Assuming that signal S ∈ R N is sparse under the sparse basis Ψ = {ψ 1 , ψ 2 · · · ψ N }, the sparse representation is where N is the length of S. The matrix form is where α is the sparse representation of S. T is a measurement matrix of size N 1 x N (N 1 N). Sampling (measurement) is a process of linear projection. The compressed signal is where Θ is a sensing matrix. It is a problem of an ill-conditioned underdetermined equation. If T and Ψ satisfy the restricted isometry property (RIP), i.e., for any constant δ k ∈ (0, 1), then K coefficients (the sparsity of S is K (K N 1 )) can be accurately reconstructed. The minimum value of δ k is the RIP constant. The equivalent of RIP is that T and Ψ are orthogonal. Signal reconstruction is an optimal solution problem of l 0 , but it is NP-hard [25]. Therefore, it is often replaced by the optimal solution problem of l 1 : where Rs is the reconstructed signal and · 1 is the l 1 norm of S. Most reconstruction algorithms are iterative greedy, convex optimization, or combination algorithms. Iterative greedy algorithms obtain the sparse vector support set of the source signal, and use the least squares estimation algorithm with limited support to reconstruct the signal [35]. Representative algorithms include matching pursuit [36] and orthogonal matching pursuit (OMP) [37]. Convex optimization algorithms convert a non-convex problem to a convex problem that is solved to approximate the source signal. Representative algorithms include basis pursuit [38] and gradient projection [39]. Combined algorithms have a high reconstruction speed, but they obtain a high probability rather than an accurate reconstruction.

Proposed Technique
In deep space exploration, our focus in the image compression of the optical sensor system is the compression (sampling) process. CS uses sparsity to perform linear compression at encoding, and nonlinear reconstruction at decoding. It effectively saves the cost of encoding and transfers the calculation cost to decoding. Therefore, it is highly suitable for remotely sensed astronomical images with a large amount of data and high redundancy in deep space exploration.
The encoding of CS has two main parts: The sparseness and the measurement matrix. The sparseness of DWT is better than that of other sparse bases, but conventional sparse vectors do not define a structure for wavelet coefficients. A method was proposed to average the wavelet coefficient matrix to construct sparse vectors [40]. Such methods do not consider the dependency between sub-bands and the direction information of highfrequency sub-bands, which will lead to poor performance in image compression. However, the fixed and improved measurement matrix cannot optimally allocate the measurement rate, and no measurement matrix is available to consider both sparse vectors and their elements. Therefore, we combine a new structured manner sparse vector and an optimized measurement matrix to design an optimal sampling scheme. The details of our proposed algorithm are as follows.

A New Sparse Vector Based on the Rearrangement of Wavelet Coefficients
The real signal in nature has sparseness in the transform domain. We can perform compression sampling for the sparse signal. The sparse representation is the essential premise and theoretical basis of CS. In this section, we construct a new sparse vector by rearranging the wavelet coefficients according to the energy distribution of the image after DWT.
For a given m × n remotely sensed astronomical image I p , the j-level DWT is applied to obtain the transformed image W I j p , which has 3j + 1 sub-bands. The size of I p satisfies m = 2 l 1 and n = 2 l 2 , where l 1 and l 2 are natural numbers. The dimension of the i-level sub-band is m 2 i × n 2 i . W I j p includes the low-frequency sub-band LL j , which reflects the approximate information of the image, and high-frequency sub-bands LH i , HL i , and HH i , i = 1, 2, · · · , j, reflecting the detailed information. Figure 2a shows the tree structure of the three-level DWT. Low-and high-frequency sub-bands LL 3 and HH 1 are at the upper-left and bottom-right corners, respectively. There is a parent-child relationship between sub-bands. The direction of the arrow is from parent to child. Except for LL 3 , HL 1 , LH 1 , and HH 1 , one parent node (m o , n o ) corresponds to four child nodes: (2m o , 2n o ), (2m o , 2n o + 1), (2m o + 1, 2n o ), and (2m o + 1, 2n o + 1). In low-frequency sub-band LL 3 , each node has only three child nodes: HL 3 , LH 3 , and HH 3 . When HL 3 , LH 3 , and HH 3 are regarded as parent sub-bands, the corresponding child sub-bands are HL 2 , LH 2 , and HH 2 , respectively. Similarly, HL 1 , LH 1 , and HH 1 are child Remote Sens. 2021, 13, 288 5 of 16 sub-bands of HL 2 , LH 2 , and HH 2 , respectively. High-frequency sub-bands HL 1 , LH 1 , and HH 1 have no child sub-band. The idea of this parent-child relationship is that the parent sub-band node and the corresponding child sub-band nodes reflect the same area information of the image. The parent sub-band has a higher importance than the child sub-band. At the same level of sub-bands, the low-frequency has the highest importance, and the importance of high frequency sub-bands, in descending order, is HL j , LH j , and HH j [41]. The parent-child relationship is the basis for calculation of the texture detailed complexity and design of the measurement matrix.
, and . When , , and are regarded as parent sub-bands, the corresponding child sub-bands are , , and , respectively. Similarly, , , and are child sub-bands of , , and , respectively. High-frequency sub-bands , , and have no child sub-band. The idea of this parent-child relationship is that the parent sub-band node and the corresponding child sub-band nodes reflect the same area information of the image. The parent sub-band has a higher importance than the child sub-band. At the same level of sub-bands, the low-frequency has the highest importance, and the importance of high frequency sub-bands, in descending order, is , , and [41]. The parent-child relationship is the basis for calculation of the texture detailed complexity and design of the measurement matrix.
We design a data structure to represent sparse vector, , according to the parentchild relationship. Figure 2b shows the first column, , of the sparse vector, , and its arrangement order is as follows (j = 3): When we design a sparse vector, we must specify the block-scanning mode. This is based on the direction information of the high-frequency sub-bands. The horizontal z-scan mode is used for , and reflects horizontal direction information; the vertical z-scan mode is used for , and reflects vertical direction information; sub-band , which reflects diagonal direction information, adopts the diagonal z-scan mode. Conventional CS sampling (measurement) is performed on each column (or each column of the sub-bands) of the coefficient matrix obtained by the DWT of the image [42,43]. High-frequency sub-bands contain the direction information of the image [44]. The performance of CS can be further improved if this characteristic can be used. Therefore, we use different scanning modes for blocks in high-frequency sub-bands. The horizontal, vertical, and diagonal z-scan modes are shown in Figure 3. The solid line in Figure 3c represents the actual scan order. We design a data structure to represent sparse vector, F, according to the parentchild relationship. Figure 2b shows the first column, F 1 , of the sparse vector, F, and its arrangement order is as follows (j = 3): where hz(·), vz(·), and dz(·) represent the horizontal, vertical, and diagonal z-scan modes, respectively. For example, dz(HH 1 ) represents the diagonal z-scan mode to arrange the coefficients of the block in HH 1 . When we design a sparse vector, we must specify the blockscanning mode. This is based on the direction information of the high-frequency sub-bands. The horizontal z-scan mode is used for HL, and reflects horizontal direction information; the vertical z-scan mode is used for LH, and reflects vertical direction information; subband HH, which reflects diagonal direction information, adopts the diagonal z-scan mode. Conventional CS sampling (measurement) is performed on each column (or each column of the sub-bands) of the coefficient matrix obtained by the DWT of the image [42,43]. Highfrequency sub-bands contain the direction information of the image [44]. The performance of CS can be further improved if this characteristic can be used. Therefore, we use different scanning modes for blocks in high-frequency sub-bands. The horizontal, vertical, and diagonal z-scan modes are shown in Figure 3. The solid line in Figure 3c represents the actual scan order. Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 15 In the new sparse vector, the number of elements contained in ( = 1,2, ⋯ , ) is where = 4 ⁄ . The new sparse vector is Table 1 shows the number of sparse vectors and their lengths for different wavelet decomposition levels.

Measurement Matrix with Double Allocation Strategy
The measurement matrix is used to observe the high-dimensional original signal to obtain the low-dimensional compressed signal. That is, the original signal is projected onto the measurement matrix to obtain the compressed signal. In this section, depending on the parent-child relationship between sub-bands and scanning modes for high-frequency sub-bands, we design a new sparse vector to achieve an optimal measurement rate allocation scheme. The new sparse vector has some unique characteristics. Each sparse vector is composed of a low-frequency sub-band coefficient and multiple high-frequency sub-band coefficients. Elements in the same sparse vector have different importance degrees. The elements of each sparse vector are the frequency information of the same area in the image. The sparse vectors have different importance to image reconstruction; image areas with rich texture are more important.
Therefore, we propose an optimized measurement matrix with a double allocation of measurement rates, which has the advantage that the effects of sparse vectors and their elements on the measurement rate allocation are considered simultaneously. For different sparse vectors, we allocate higher measurement rates to important sparse vectors. In a sparse vector, more weight is given to important elements. The main steps are as follows.
1. An × 4 ( < 4 ) random Gaussian matrix is constructed. 2. In a sparse vector, the importance of the elements decreases from the front to the back. Increasing the coefficient of the first half of the measurement matrix can preserve more important information of the image. To do this, and to satisfy the incoherence requirement of any two columns of the measurement matrix, we introduce an × optimized matrix, In the new sparse vector, the number of elements contained in F i (i = 1, 2, · · · , l) is where l = mn/4 j . The new sparse vector is Table 1 shows the number of sparse vectors and their lengths for different wavelet decomposition levels.

Measurement Matrix with Double Allocation Strategy
The measurement matrix is used to observe the high-dimensional original signal to obtain the low-dimensional compressed signal. That is, the original signal is projected onto the measurement matrix to obtain the compressed signal. In this section, depending on the parent-child relationship between sub-bands and scanning modes for high-frequency subbands, we design a new sparse vector to achieve an optimal measurement rate allocation scheme. The new sparse vector has some unique characteristics. Each sparse vector is composed of a low-frequency sub-band coefficient and multiple high-frequency subband coefficients. Elements in the same sparse vector have different importance degrees. The elements of each sparse vector are the frequency information of the same area in the image. The sparse vectors have different importance to image reconstruction; image areas with rich texture are more important.
Therefore, we propose an optimized measurement matrix with a double allocation of measurement rates, which has the advantage that the effects of sparse vectors and their elements on the measurement rate allocation are considered simultaneously. For different sparse vectors, we allocate higher measurement rates to important sparse vectors. In a sparse vector, more weight is given to important elements. The main steps are as follows. 1. An In a sparse vector, the importance of the elements decreases from the front to the back. Increasing the coefficient of the first half of the measurement matrix can preserve more important information of the image. To do this, and to satisfy the incoherence where z 1 and z 2 are different prime numbers. Since z 1 and z 2 are relatively prime, any two columns of W i are incoherent [45,46].

3.
The upper left m i × m i part of G i is dot-multiplied with W i , and other coefficients of G i remain unchanged, to obtain an optimized measurement matrix G i . Because z 1 > 1 and z 2 > 1, the coefficients of the upper-left m i × m i part of G i are increased. 4. m i is the size of W i . It is determined by the measurement rate of each sparse vector. When the total measurement rate is certain, the measurement rate is allocated according to the texture detail complexity of the image. We calculate m i as follows.

•
The elements in the sparse vector correspond to the frequency information of the same area in the image. The high-frequency coefficients reflect the detailed information. The energy of high-frequency coefficients is defined to describe the texture detailed complexity of the sparse vector. Taking F 1 as an example, the energy of the high-frequency coefficients of F 1 is where F 1 (i) represents the i-th element in F 1 .

•
If the total measurement rate is R total , then the total measurement number is N total is divided into two parts: Adaptive distribution and fixed distribution. The measurement number used to adaptively allocate according to the complexity of the texture is The measurement number of fixed allocation is i.e., N 2 is evenly allocated to each sparse vector.

•
The higher the texture detailed complexity, the higher measurement rates are allocated to ensure retention of more details. Image areas with low texture detailed complexity, such as the background and smooth surface areas, are allocated lower measurement rates. The measurement rate is adaptively allocated according to the texture detailed complexity. A linear measurement rate allocation scheme is established based on this principle. The measurement number of the sparse vector F i is where mn/4 j is the number of sparse vectors, and is a proportionality coefficient. Then the measurement number of F i is and the final measurement number of F i is round(NF i ), i.e.,

5.
The sparse vectors are compressed by the optimized measurement matrix, and the compressed value of F i is

Image Reconstruction
The reconstruction is an optimal solution problem for solving l 1 , where RF i is the recovered sparse matrix, s.t. is the abbreviation of "subject to". We use the OMP algorithm to reconstruct sparse vectors [25]. OMP is a classical greedy iterative algorithm that follows the atomic selection strategy of the matching tracking algorithm. Each iteration selects an atom, and the selected atomic set should be normalized in each iteration to ensure the optimization of each iteration result, thus reducing the number of iterations and accelerating convergence. OMP performs vector orthogonalization during each iteration. OMP has higher space and time complexity than the matching pursuits method. In deep space exploration, our optical sensing system compresses remotely sensed astronomical images and transmits the compressed data, and associated calculation costs, to the ground. The recovered vectors are rearranged into a discrete wavelet coefficient matrix, and the reconstructed image is obtained through an inverse discrete wavelet transform.

Evaluation Standard
The peak signal to noise ratio (PSNR) and structural similarity (SSIM) [47] are employed as evaluation indexes of the reconstructed image quality. The PSNR is widely used to measure the distortion effect after compression. As a quantitative measure to the subjective quality, the SSIM demonstrates the quality of detailed information such as the textures and edges of the reconstructed image [48].
The PSNR of an m × n image is where is the mean square error, and I and RI, respectively, represent the original and reconstructed image. For a color image, where L = 1, 2, 3 stands for the three channels in the image. In this paper, we converted RGB images to YUV channels for processing. The SSIM evaluates the luminance, contrast, and structure of the image. It can be represented as: where µ I , µ RI , σ 2 I , σ 2 RI , and σ IRI are the local mean, variance, and cross-covariance for I and RI.
where Y is the maximum value of pixel gray, δ 1 = 0.01, δ 2 = 0.03. For a color image, the SSIM value is the average of the three channels.

Experimental Data and Evaluation
Experiments were performed on a simulated celestial body image, the gibbous moon, and images of planets, and the experimental results were evaluated subjectively and objectively by comparing three algorithms. Image compression in deep space detection is mainly based on JPEG2K, an algorithm with wavelet transform as the core technology. Starck et al. [29] analyzed the application and advantages of traditional CS algorithms on astronomical images. Traditional CS algorithms compress each column of the sparse matrix after DWT and use a random Gaussian measurement matrix with a fixed measurement rate. A DWT-CS algorithm proposes [40] to construct a sparse vector, giving more weight to higher energy coefficients. Three-level DWT is performed on an image using a coif5 wavelet. The prime numbers z 1 and z 2 of the measurement matrix are 2 and 3, respectively. The images of moon and planets are taken from Celestial, which is an open-source software for astronomical works supported by NASA Technology [49]. The PSNR and SSIM are the average of 1000 experiments.

Simulated Images of Celestial Bodies
A simulated celestial body image can be obtained by calculating the number of photoelectrons on each illuminated pixel [10]. The number of photoelectrons is where Q e is the quantum efficiency of the image sensor, F 0 is the fill factor, S pixel is a single pixel area, ∆t is the exposure time, and the total photon flux hitting the pixel is where R l is the radiance of the reflected light, τ is the transmittance of the optics, A d is the aperture diameter, f 0 is the focal length, and θ 0 is the angle between the line-of-sight vector and the observer surface normal [10].
The celestial body imaging model is the convolution of the cumulative energy distribution of the optical sensor system with the PSF [3], where PSF(x, y) is the point spread function of the optical system, and * is the convolution operation.
We simulated a 256 × 256 celestial body image whose radius was 100 pixels. where ( , ) is the point spread function of the optical system, and * is the convolution operation. We simulated a 256 × 256 celestial body image whose radius was 100 pixels. Figure 4 displays its normalized intensity distribution model. Tables 2 and 3, respectively, show the values of PSNR and SSIM when CR is 0.1, 0.2, 0.3, 0.4, 0.5, and 0.6. In this paper, the CR is the ratio of the compressed image to the original image.   When CR is low, the CS, DWT-CS, and the proposed algorithm perform better than JPEG2K because CS takes advantage of sparsity and can restore the original data through a lower sampling rate. As CR increases, the reconstruction quality of JPEG2K improves significantly. JPEG2K achieves low CR and high fidelity with difficulty, resulting in a long transmission time or a poor image quality [50]. The proposed algorithm performs best when CR is 0.1-0.6. It constructs a new sparse vector based on the energy distribution of the image in the frequency domain and the information carried by the sub-bands, and double allocation of the measurement rate can realize an optimal sampling scheme. When the total measurement rate is determined, more measurement weight is given to the sparse vector with richer surface texture and elements with higher energy, ensuring better reconstruction quality at a lower CR.

Moon Image
As the starting point of human deep space exploration, the moon has great significance. Figure 5 shows the moon image and a partial enlargement. The moon has abundant features, such as texture, edges, and lunar mare.

Moon Image
As the starting point of human deep space exploration, the moon has great cance. Figure 5 shows the moon image and a partial enlargement. The moon has ab features, such as texture, edges, and lunar mare.  Table 4 compares the PSNR values of the four algorithms. The CR ranges fro 0.6, and the interval is 0.1. The PSNR values of the four algorithms increase with CR pared to JPEG2K, the PSNRs of the proposed algorithm improve by 7.2812 and 1.   Table 5 shows the SSIM values of the four algorithms. When CR is 0.1, compared with JPEG2K, CS, and DWT-CS, the SSIM values of the proposed algorithm are improved by 0.1351, 0.0877, and 0.0583, respectively. When CR is 0.1-0.6, the SSIM of the proposed algorithm show the best performance, which indicates that our algorithm can retain more texture information at the same CR compared with other algorithms.  Figure 6 shows part of the reconstructed image when the CR is 0.1 and 0.6. When the CR is low, the reconstructed image of JPEG2K has obvious distortion. With the increase of CR, the quality of the reconstructed images by the four algorithms improves significantly. Comparing the edge and surface textures of the reconstructed image, the proposed algorithm performs better than the others. The scanning mode based on direction information and the measurement rate allocation strategy shows a significantly improved reconstruction effect.
CR is low, the reconstructed image of JPEG2K has obvious distortion. With the increase of CR, the quality of the reconstructed images by the four algorithms improves significantly. Comparing the edge and surface textures of the reconstructed image, the proposed algorithm performs better than the others. The scanning mode based on direction information and the measurement rate allocation strategy shows a significantly improved reconstruction effect.  We analyze the average PSNR and SSIM of the three planets. Figure 8 shows PSNR values of the three planets when CR is 0.1-0.6 (with an interval of 0.1), and all crease with CR. The proposed algorithm has evident advantages when CR is 0.1, w PSNR values increasing by 6.3291, 4.2989, and 2.9640 dB, respectively, compared JPEG2K, CS, and DWT-CS. The values of the algorithms are relatively close when CR 0.6, increasing by 2.8765, 2.7942, and 1.6371 dB for the proposed algorithm over JPEG2 CS, and DWT-CS, respectively. Table 6 shows the SSIM values of the three planets. Wh CR is 0.1-0.6, compared with the other three algorithms, the SSIM of the proposed al rithm have better performance. When CR is 0.1, compared with JPEG2K, CS, and DW CS, the SSIM values of the proposed algorithm are improved by 0.1010, 0.0839, and 0.05 respectively. When CR is low, the reasonable use of the measurement rate is particula important. Our double allocation scheme can optimally allocate the limited measurem rate to obtain a better reconstructed image.   We analyze the average PSNR and SSIM of the three planets. Figure 8 shows the PSNR values of the three planets when CR is 0.1-0.6 (with an interval of 0.1), and all increase with CR. The proposed algorithm has evident advantages when CR is 0.1, with PSNR values increasing by 6.3291, 4.2989, and 2.9640 dB, respectively, compared to JPEG2K, CS, and DWT-CS. The values of the algorithms are relatively close when CR is 0.6, increasing by 2.8765, 2.7942, and 1.6371 dB for the proposed algorithm over JPEG2K, CS, and DWT-CS, respectively. Table 6 shows the SSIM values of the three planets. When CR is 0.1-0.6, compared with the other three algorithms, the SSIM of the proposed algorithm have better performance. When CR is 0.1, compared with JPEG2K, CS, and DWT-CS, the SSIM values of the proposed algorithm are improved by 0.1010, 0.0839, and 0.0521, respectively. When CR is low, the reasonable use of the measurement rate is particularly important. Our double allocation scheme can optimally allocate the limited measurement rate to obtain a better reconstructed image. Figure 9 shows a partially enlarged view of the reconstructed images when CR is 0.1. We analyze the edges of Jupiter and Mars. The image edges reconstructed by the four algorithms are all distorted. In particular, for JPEG2K, the edge distortion is more serious, and the blocking effect is evident. Our algorithm also has distortions that can be distinguished, but it performs better than the other algorithms. Impact craters are used to analyze the reconstruction effect of Mercury. Compared to the original image, the other algorithms have obvious distortions, and the reconstruction effect of our algorithm is closest to the original image. Our algorithm preserves the shape of the impact crater and obtains better reconstruction quality.
CS, and DWT-CS, respectively. Table 6 shows the SSIM values of the three planets. When CR is 0.1-0.6, compared with the other three algorithms, the SSIM of the proposed algorithm have better performance. When CR is 0.1, compared with JPEG2K, CS, and DWT-CS, the SSIM values of the proposed algorithm are improved by 0.1010, 0.0839, and 0.0521, respectively. When CR is low, the reasonable use of the measurement rate is particularly important. Our double allocation scheme can optimally allocate the limited measurement rate to obtain a better reconstructed image.     Figure 9 shows a partially enlarged view of the reconstructed images when CR is 0.1. We analyze the edges of Jupiter and Mars. The image edges reconstructed by the four algorithms are all distorted. In particular, for JPEG2K, the edge distortion is more serious, and the blocking effect is evident. Our algorithm also has distortions that can be distinguished, but it performs better than the other algorithms. Impact craters are used to analyze the reconstruction effect of Mercury. Compared to the original image, the other algorithms have obvious distortions, and the reconstruction effect of our algorithm is closest to the original image. Our algorithm preserves the shape of the impact crater and obtains better reconstruction quality.

Conclusions
We presented a wavelet-based sensing algorithm to compress remotely sensed astronomical images, fully considering the characteristics of an image in the frequency domain and using a sparse matrix construction method based on the parent-child relationship between sub-bands and different scanning modes of high-frequency sub-bands. The optimized measurement matrix with a double allocation of measurement rates provides an

Conclusions
We presented a wavelet-based sensing algorithm to compress remotely sensed astronomical images, fully considering the characteristics of an image in the frequency domain and using a sparse matrix construction method based on the parent-child relationship between sub-bands and different scanning modes of high-frequency sub-bands. The optimized measurement matrix with a double allocation of measurement rates provides an optimal measurement rate allocation strategy. The optimized measurement matrix retains the most important information of the frequency domain at a low measurement rate. An astronomical image was reconstructed by OMP and IDWT. The proposed algorithm was verified on remotely sensed astronomical images. Experiments demonstrated its advantages in PSNR and SSIM compared to JPEG2K, CS, and DWT-CS. We have provided a high-performance remotely sensed astronomical image compression algorithm for a miniaturized independent optical sensor system in deep space exploration. The algorithm has a certain reference significance for other images with a rich texture.
Author Contributions: Y.Z. is responsible for the research ideas, overall work, the experiments, and writing of this paper. J.J. provided guidance and modified the paper. G.Z. is the research group leader who provided general guidance during the research and approved this paper. All authors have read and agreed to the published version of the manuscript.