FUSION OF MULTI-SCALE DEMS FROM DESCENT AND NAVCM IMAGES OF CHANG ’ E-3 USING COMPRESSED SENSING METHOD

The multi-source DEMs generated using the images acquired in the descent and landing phase and after landing contain supplementary information, and this makes it possible and beneficial to produce a higher-quality DEM through fusing the multi-scale DEMs. The proposed fusion method consists of three steps. First, source DEMs are split into small DEM patches, then the DEM patches are classified into a few groups by local density peaks clustering. Next, the grouped DEM patches are used for subdictionary learning by stochastic coordinate coding. The trained sub-dictionaries are combined into a dictionary for sparse representation. Finally, the simultaneous orthogonal matching pursuit (SOMP) algorithm is used to achieve sparse representation. We use the real DEMs generated from Chang’e-3 descent images and navigation camera (Navcam) stereo images to validate the proposed method. Through the experiments, we have reconstructed a seamless DEM with the highest resolution and the largest spatial coverage among the input data. The experimental results demonstrated the feasibility of the proposed method. * Corresponding author.


INTRODUCTION
Chang'e-5 will be China's first lunar sample return mission.High precision topographic mapping of Chang'e-5 landing site can provide detailed terrain information to ensure the safety of the lander as well as to support tele-operated sampling of lunar soils and rocks.The lander will acquire descent images in the descent and landing phase, and will also acquire stereo images of the sampling area after landing.A variety of Digital Elevation Model (DEM) products will be generated using image data acquired by different sensors in different phases and they often have differences in spatial coverage, resolution, and accuracy.The multi-source DEMs contain supplementary information, and this makes it possible and beneficial to produce a higher-quality DEM through fusing the multi-scale DEMs.
In the previous Mars and lunar landed missions, descent images have been used for localization of the lander and/or mapping of the landing area.Li et al. (2002) conducted experiments with simulated descent images and Field Integrated Design and Operations Rover data.Ma et al. (2001) developed an integrated bundle adjustment system incorporating both descent and roverbased images to localize the rover along the traverse.Approaches to visual localization using overhead and ground images were reviewed and the particular capabilities under development at JPL were discussed (Matthies et al., 1997).In Mars Exploration Rover (MER) mission, three descent images of low resolution were acquired by Spirit and Opportunity respectively, and they were used in lander localization but was not further used for 3D mapping of the landing area (Li et al., 2005).In the Mars Science Laboratory (MSL) mission, the descent images acquired by the Mars Descent Imager (MARDI) in the Entry-Descent-Landing (EDL) phase were compiled into image mosaics and provide color coverage of the landing site and science target regions; the mosaics are incorporated into the landing base map (Parker et al., 2013).In Chang'e-3 (CE-3) mission, descent images were used to generate high precision topographic products of the landing site with different resolutions (0.05 m, 0.2 m,0.4 m) and was also used in lander localization (Liu et al., 2015).
For planetary exploration as well as earth observation applications, high spatial resolution DEMs are always of limited spatial coverage due to the high cost of data acquisition, and may have data quality problems such as data voids and noises, while relatively low-resolution DEMs provide less spatial information but usually cover larger areas.In the past years, many studies have been conducted on improving the quality of DEMs, among which fusion ideas haven been introduced into DEM reconstruction.For example, SRTM and ASTER GDEM data were fused in the frequency domain, and the data voids were filled, so that the overall accuracy of the fused data was improved (Karkee et al., 2008).Multi-scale modeling is adopted to fill the voids in high resolution DEM data, and a multi-scale Kalman smoother (MKS) was used to remove blocky artifacts in DEM fusion (Jhee et al., 2013).These methods are not applicable for fusing DEMs with different resolutions, coverage, and vertical accuracies.
To address inhomogenity of available DEM products, several methods of fusing DEMs have been developed to obtain a complete DEM coverage with improved quality.Recently, sparse-representation based methods, as a subset of transformdomain fusion methods, have been applied to DEM fusion.Papasaika et al. (2011) presented a generic algorithmic approach for fusing two arbitrary DEMs, using the framework of sparse representations, and experiments with real DEMs from different earth observation satellites were conducted.Yue et al. (2015) proposed a regularized framework for production of high resolution DEM data with extended coverage.Boufounos et al. (2011) introduced a new sparsity model for fusion frames, and the model provides a promising new set of mathematical tools and signal models useful in a variety of applications, probabilistic analysis shows that under very mild conditions the probability of recovery failure decays exponentially with increasing dimension of the subspaces.Divekar and Ersoy (2009) created a dictionary that relates high resolution image patches from a panchromatic image to the corresponding filtered low resolution versions, and proposed two algorithms which directly use the dictionary and its low resolution version to construct the fused image.Tao and Qin (2011) proposed a new image fusion algorithm in the compressive domain by using an improved sampling pattern.One key advantage offered by the compressed sensing (CS) approach is that samples can be collected without assuming any prior information about the signal being observed, thereby motivating our research on compressive fusion of DEM.
In this paper we present an improved CS method for fusing multi-scale DEMs to produce high-resolution DEM with extended coverage.Through the experiment using Chang'e DEMs, we have reconstructed a seamless DEM with the highest resolution and the largest spatial coverage among the input data.
The rest of this paper is structured as follows: Section 2 provides a brief description of compressed sensing; Section 3 presents and specifies the proposed method; Experimental results are presented in Section 4. Finally, conclusions and suggestions for future work are given in Section 5.

COMPRESSED SENSING
Compressed sensing is a signal processing technique for efficiently acquiring and reconstructing a signal, by finding solutions to underdetermined linear systems.This is based on the principle that, through optimization, the sparsity of a signal can be exploited to recover it from far fewer samples than required by the Shannon-Nyquist sampling theorem (Donoho, 2006).
Supposing the signal f can be recovered from a set of M measurements.This compressive measurement vector can be formulated as where is a measurement matrix.Since M<<N, the recovery of the signal vector f from the measurement vector y is a highly underdetermined problem.However, there are two conditions under which recovery is possible.The first one is sparsity which requires the signal to be sparse in some domain.The signal f can be represented sparsely by an orthogonal basis.The second one is incoherence which is applied through the isometric property which is sufficient for sparse signals.The orthogonal basis  and compressive measurement matrix are incoherent.

METHOD
The proposed fusion method consists of three steps.First, source DEMs are split into small DEM patches, then the DEM patches are classified into a few groups by local density peaks clustering.Next, the grouped DEM patches are used fordictionary learning by K-SVD algorithm.Finally, the simultaneous orthogonal matching pursuit (SOMP) algorithm is used to achieve sparse representation.After the three steps, the resultant sparse coefficients are then fused following the max L1-norm rule.The fused coefficients can be inversely transformed to a high-resolution DEM with extended coverage by using the learned dictionary.

Problem formulation
Supposing yk representing DEM of different coverage and resolution, the generative model of multi-DEM can be defined as where h  represents noise vectors; Lk is set as the downsampling and operator; Mk represents the transformation matrix; as the coverage for each DEM differs, the cropping operator Ok is defined as a diagonal matrix with zero elements if the pixel was invalid in the kth input data; furthermore, the voids and the anomalies in the DEM are also included in Ok.The problem is to fuse measurements yk(k=1,2,…n) to recover DEM x.
Assuming the fused result x can be represented as sparse linear combination of elements from a dictionary D, while the elements of the D are called atoms.The result x is sparsely represented over D if x=Da0, a0 denotes a sparse coefficient vector with most zeros.Therefore, the generative model can be represented as where Dk is defined as dictionaries of different resolution.
Given Dk and yk are available, the sparse coefficients a0 can be recovered, and the fused DEM x can determined by computing 0 k D  .The function can be expressed as The first term corresponds to the reconstruction error with respect to the observed DEMs yk.The second term is associated with the L1 norm of the candidate solution vector α.The parameter τ controls the trade-off between data fitting and sparsity.However, since each point at different DEMs has different accuracy, weights parameters should be included in the problem formulation.Therefore the optimization function is modified as follows: where wk denotes the weight for the kth DEM.
In order to deal with the blocking artifacts along the patch borders, the consistency between neighboring patches is imposed.Assuming operator P extracts the overlap region between patches, yp is a vector containing DEM values in the overlap region, adding the regularization term that minimizing the discrepancy between overlapping patches into (4), the final formulation can be expressed as where parameter β control the influence of the patch overlap factor.

Dictionary construction
As learning an over-complete dictionary capable of representing classified of DEM patches is of great importance.The raw patches are randomly sampled in this paper, which is similar to the approach used in Yang et al.(2008).Then the sampled DEM patches were classified into several new groups based on feature vectors.Since K-SVD is one of the most popular dictionary learning algorithms, the sub-dictionaries are trained by K-SVD.
As the structures of DEM patches in each cluster are similar, sub-dictionary learning scheme can get more accurate structure description of input DEM patches.The detail of K-SVD based sub-dictionary learning and combination are shown as follows: Step 1: A few sub-dictionaries S1, S2, ..., Sn are learned for each DEM patch groups Step 2: Sub-dictionaries of each cluster are trained by K-SVD algorithm.Then the sub-dictionaries are combined to a dictionary for fusion.where  is the combined dictionary and S1, S2, …, Sn are the trained sub-dictionaries.

Adaptive-weight determination for DEM fusion
The fusion is completed with weight maps that reflect the estimated relative accuracy of the source DEMs at each grid point.First, geometric registration should be for the datasets.The transformation matrix Mk can be obtained after registration, and the cropping region can also be derived according to the coordinates.
A data-driven strategy is used to find the weights based on geomorphological characteristics, which is similar to the approach used Boufounos et al. (2011).Based on the slope and entropy, the resulting accuracy maps of both input DEMs at each overlapping point can be determined, finally the reciprocal values are used as weights for the fusion (Papasaika et al., 2011).

Fusion of DEMs
Since optimization problems of this form constitute the main computational kernel of compressed sensing applications, there exists a wide options of algorithms for their solution.As its simplicity and computational efficiency, Orthogonal Matching Pursuit (OMP) (Mallat,1998) is adopted.In the experiment, the overlap parameter β is set between 0.6 and 1.2.The number of non-zero atoms is set between 7 and 15.The minimum patch size is set as 3×3 and it should not be bigger than 9×9.

Landing site mapping using CE-3 descent images
CE-3 began to descent from the lunar orbit at an altitude of around 15 km, and when it was about 2 km above the lunar surface, the descent camera started to take images.During the phases of descending, hovering and obstacle avoidance and landing, the descent camera acquired totally 4,672 images with a resolution higher than 1 m within an area of 1*1 km and as high as 0.1 m within a range of 50 m from the landing point.Main technical parameters of the CE-3 descent camera are listed in table 1.One hundred and eighty were selected with equal time interval and incorporated in a self-calibration free network bundle adjustment, and the initial trajectory (including camera positions and attitudes) of the camera was recovered.Then, 26 ground control points (GCPs) are selected from the rectified Chang'e-2 DEM and digital orthophoto map (DOM) for absolute orientation.The RMSEs (Root Mean Square Errors) of these GCPs are 0.724 m, 0.717 m and 0.602 m in three directions.The RMSEs of 18 check points are also less than 1 m. Figure 1 shows the DEM and DOM generated from 80 descent images with a resolution of 0.075 m.The dots in the maps represent the lander position.The maps cover an area of 97 m  115 m.Using more descent images of higher altitudes, topographic products of larger coverage were also generated.

Image size
Actual imaging distance

Mapping with Navcam images
The geometric parameters of Yutu rover's Navcam are listed in Table 2. Figure 2 shows DEM at waypoint D, it has a resolution of 0.02 and was automatically generated from the Navcam images.In order to fuse DEMs with descent images, the ground DEM is transformed into the lunar body-fixed coordinate system.

Integration of mapping products
Discrepancies exist between the topographic products from descent and ground images, co-registration between the two DEMs were performed using Integrative Closet Point algorithm.
After the co-registration, the mean differences between the two data sets have been reduced from 0.142 m, 0.096 m, 0.765 m to 0.003 m, 0.004 m, 0.221 m in three directions respectively, as shown in Figure 3.
As in Equation ( 8), K represents the reconstructed measurement, while I is the reference data, I MAX is the maximum possible pixel value of the image.Smaller PMSE and larger PNSR correspond to better performance.Figure 6 shows that our proposed method has better performance and can reconstruct an enhanced elevation result.The quantitative results in

CONCLUSIONS
In this paper we present an improved compress sensing method for fusing multi-scale DEMs to produce high-resolution DEM with extended coverage.Compared with traditional fusion methods, the reconstructed data are generated using the supplementary information between DEMs with different resolution and coverage.Furthermore, the grouped DEM patches instead of the whole source image are used fordictionary learning in the proposed method, which can improve efficiency of the method.We use the real DEMs generated from Chang'e-3 descent images and Navcam stereo images to validate the proposed method.Through the experiments, we have reconstructed a seamless DEM with the highest resolution and the largest spatial coverage among the input data.The experimental results demonstrated the feasibility of the proposed method.However, there were still some limitations to the proposed method.In the future, more complementary factors like the edginess, the noisiness of DEM will be taken into account to improve the accuracy of the fused data.
Figure 1.DEM (right) and DOM (left) generated from of descent images

Figure 3 .
Figure 3. Co-registration of two datasets

Figure 6 .Figure 8 .
Figure 6.The reconstruction results of bilinear, kriging interpolation, and the proposed method.

Table 3
also confirm the trend.

Table 3 .
The quantitative results