Super-Resolution Geomagnetic Reference Map Reconstruction Based on Dictionary Learning and Sparse Representation

Geomagnetic matching navigation plays a vital role in the navigation and guidance field, and the construction of geomagnetic reference maps plays a significant role in geomagnetic matching navigation. This paper addresses the problem of generating high-resolution geomagnetic reference maps with limited measured data. Most existing methods can hardly satisfy the accuracy requirements. We base our study on the theory of image super-resolution reconstruction and approach this problem from the perspective of dictionary learning and sparse representation. We propose a method of sparse dictionary initialization based on prior information from rectangular harmonic analysis, and then the K-SVD algorithm is applied to the dictionary training process to improve the performance of the sparse dictionary. Three components of the geomagnetic field are considered for multi-channel sparse representation to enhance the quality of the constructed maps. Experimental results show that our proposed method outperforms other geomagnetic reference map construction methods in the reconstructed precision as well as the robustness to noise.


I. INTRODUCTION
The mainstream development direction for aircraft navigation technology in recent years has been the satellite navigation and inertial navigation technology [1]- [3]. However, the anti-interference ability of satellite navigation is poor, and the error of inertial navigation accumulates over time. The geomagnetic field changes slowly with time and is not affected by climate and region. Hence, geomagnetic navigation technology has gradually attracted attention as an auxiliary navigation method owing to its strong anti-interference ability [4], [5].
Geomagnetic matching navigation technology uses realtime acquired geomagnetic field data to match the preprepared geomagnetic reference map of the target region to determine the location of vehicles [6]. Therefore, the accuracy of the geomagnetic reference map construction plays a decisive role in the performance of geomagnetic navigation [7]. The steps of geomagnetic navigation are shown in Fig. 1. Firstly, the geomagnetic information of target areas is The associate editor coordinating the review of this manuscript and approving it for publication was Michele Nappi . measured to build geomagnetic reference maps in advance. Secondly, flight vehicles measure the real-time geomagnetic strength when they fly over these areas. Finally, they can be positioned by matching their real-time geomagnetic strength with the geomagnetic maps. There are two main methods for constructing geomagnetic reference maps. One is based on the existing physical model of the geomagnetic field, and the other is constructing a gridded geomagnetic reference map based on measured geomagnetic data [8]. The existing International Geomagnetic Reference Magnetic Field (IGRF) and World Magnetic Model (WMM) describe the earth's primary magnetic field [9]. Since the geomagnetic field model relies on satellite monitoring data and the anomalous field inside the earth is attenuated as the altitude increases, the geomagnetic field model contains a significant deviation from real-time measured geomagnetic data. Geomagnetism based on measured data is more reliable than the geomagnetic field model. Previous work has only focused on interpolation methods for geomagnetic reference map construction. Commonly used methods include cubic spline interpolation, Kriging interpolation, and particle swarm optimization Kriging (PSO-Kriging) interpolation [10]. Although the existing algorithm can improve the peak signal-to-noise ratio of the reference map, it is easy to produce an excessively smooth result, and it is difficult to recover the details after smoothing. Therefore, it is essential to improve the resolution of the geomagnetic reference map.
Sparse representation has the potential to extract the textures of signals and is widely used in image processing. Inspired by prior image super-resolution via sparse, we have proposed a high-precision geomagnetic reference mapbuilding methodology based on sparse representation and dictionary learning (SR-GRM). The main contributions of this paper are as follows: (1) The design and implementation of a high-precision geomagnetic reference map-building methodology based on sparse representation and dictionary learning.
(2) The method of sparse dictionary initialization based on prior information from rectangular harmonic analysis and experimental results indicate that the proposed method can generate a more qualitative sparse dictionary with a small training set size.
(3) The precision of geomagnetic reference map reconstruction for SR-GRM with bicubic interpolation, neighbor embedding, and PSO-Kriging interpolation are compared.
(4) The noise robustness of the proposed method is compared with different kinds of other super-resolution methods. The experimental results illustrate that our method can outperform other methods in reducing the impact of noise. This paper is organized as follows. In section II, we review the relevant research in the field of single image superresolution, especially sparse representation, and dictionary learning. A detailed description of the proposed SR-GRM algorithm is discussed in section III. In section IV, experimental results demonstrate the efficacy of SR-GRM as a priori for geomagnetic reference map building. Finally, some concluding remarks are given in section V.

II. RELEVANT WORK
In the image super-resolution field, all methods can be classified into one of three categories: interpolation-based methods, reconstruction-based methods, or learning-based methods [11].
Interpolation-based methods have the advantage of reducing the amount of computation and are recognized as the most classical algorithms in image resolution enhancement. However, these methods do not contribute to generating details between pixels in low-resolution (LR) images, and additional artifacts are unavoidable.
The central concept in construction-based methods is that the same LR images can be produced by applying the image formation model to the reconstructed image. Nevertheless, much detailed information is lost when generating LR images from high-resolution (HR) images [12]. Therefore, the recovery problem is underdetermined, and each LR image may have several corresponding solutions. Various solutions have been proposed for the problem. For example, one can choose a maximum a-posteriori (MAP) solution under generic images, such as the Markov Random Field (FMR) [13]- [15]. In [16], a multi-modality image fusion method based on image decomposition and sparse representation is proposed, which shows strong representation ability to the training set.
Recently, inspired by the rapid development of deep learning in computer vision technology, researchers began to use neural networks in the single-image super-resolution field [17], [18]. Deep convolutional networks (CNN) [19], [20] and generative adversarial networks (GAN) [19] are designed to learn the non-linear relationship between the LR datasets and HR detests. However, enormous datasets containing millions of LR and HR patch pairs are required to guarantee the datasets are expressed sufficiently. Concerning the difficulty in collecting geomagnetic information, it is impractical to build satisfactory datasets for neural network training.
Existing interpolation-based methods such as PSO-Kriging interpolation and bicubic are easy for calculation but influence the precision of reconstructed maps. Learning-based methods can only show its superiority with a sufficiently big dataset for training. However, existing geomagnetic data can hardly satisfy the demand of training set. Compared with the learning-based methods described above, a much smaller dataset is required in the proposed sparse representation algorithm. Only the LR dictionary is used in the recovery process of sparse representations. As the computation of our algorithm is based on linear regression, computing efficiency is greatly improved. Besides, the K-SVD [21]- [23] algorithm is applied to the dictionary training process, which can maintain the precision of sparse representation and computational efficiency at the same time.

III. PROPOSED METHODOLOGIES
Three of the seven elements in the geomagnetic field are independent, so three dimensions need to be considered separately when building the geomagnetic reference maps [2]- [4].  In this paper, the three independent components of the geomagnetic field are chosen as geomagnetic north, geomagnetic east, and vertically downward, and are denoted as X, Y, and Z, respectively. Fig. 2 shows the features of the geomagnetic field. Since the total field strength is the main component of interest in geomagnetic navigation, this paper uses the three independent components to reconstruct high-resolution geomagnetic reference maps and then combines them into the total field strength for performance evaluation of our method. Fig. 3 shows the framework of our solution which includes the three procedures below. (1) Geomagnetic information from the United States Geological Survey (USGS) and International Geomagnetic Reference Field (IGRF) is prepared as data set. Dictionary learning method is used to extract features of the data set. (2) The LR reference maps are constructed by meshing the real-time measured data through Kriging interpolation method. (3) The HR reference maps are reconstructed based on the theory that they have the same sparse codes of given features with the corresponding LR reference maps.
Based on the discussion above, geomagnetic reference maps can be constructed using the concept of super-resolution image reconstruction. The algorithm proposed in this paper includes the following two parts: sparse representation and dictionary learning.

A. SINGLE CHANNEL SPARSE REPRESENTATION
In this paper, the low-resolution reference map, Y , is obtained by blurring and downsampling the high-resolution map X (Eq. (1)).
where H is the downsampling matrix, and L is the blurring filter. It can be seen from Eq. (1) that super-resolution reconstruction is an ill-conditioned problem [21]. Therefore, Eq.
(1) is an underdetermined system, and each map corresponds to an infinite number of solutions. Nevertheless, the problem can be solved well under the following constraints. The first constraint is that the small patches x of the highresolution map X can be sparsely represented by a redundant dictionary D h , which can be formulated according to Eq. (2). where α is the sparse representation coefficient of the HR patch x. Since the LR patch and the corresponding HR patch have the same textures, two dictionaries D l and D h need to be trained, and the LR and the corresponding HR patches can be sparsely represented with the same coefficients. In other words, a sparse representation of input LR patch y is obtained with respect to D l , and the corresponding output HR patch x will be generated by multiplying D h by α. The sparsest representation of the HR patch y can be acquired from Eq. (3).
where F is a feature transformation matrix. It is used to extract the high-frequency components of the LR map. Since the high-frequency content can better distinguish the target area from others, F works as a high-pass filter to rebuild the lost high-frequency components in the HR map. Although the l 0 -norm minimization problem in (3) is NP-hard [21], it can be effectively recovered by minimizing the l 1 -norm while α are sufficiently sparse. Equation (3) can be reformulated as Eq. (4).
where λ balances the residue of the sparse representation to y and sparsity of coefficients α. Equation (5) is essentially a convex optimization problem and has been excellently solved. Solving the problem in (5) alone does not guarantee consistency between adjacent map patches. Therefore, when the formula (5) is satisfied, it is also necessary to add a restriction condition that the overlapping portions between the reconstructed adjacent high-resolution image blocks are to be consistent. The optimized problem can be expressed by Eq. (6).
where P is the matrix for extracting the overlapping region between the previously reconstructed HR map and the current patch, and ω is composed of the overlapping part of the previously reconstructed HR map. The optimization in Eq. (6) can be reformulated as Eq. (7).
where D = FD l βPD h , and y = Fy βω . β is set to balance the proximity of the low-resolution input and compatibility with its neighbors. The parameter β is set to 1 in our experiments. An optimal high-resolution map can be reconstructed by calculating the solution α * from Eq. (7). The final optimization problem can be formulated as Eq. (8).
where X * is the optimal reconstructed high-resolution map. The proposed algorithm of sparse representation is shown in Fig. 4.

B. MULTI-CHANNEL SPARSE REPRESENTATION
Considering that the high-frequency parts of the HR patches of the three geomagnetic signal channels are similar, the correlation between the three channels may be enforced according to Eq. (9).
S µ x µ − S γ x γ 2 2 < ε µν , µ, γ ∈ {n, e, v} , µ = γ , (9) where n, e, and v indicate the three channels of the geomagnetic signal, and the S matrix is a high-pass filter. For example, S n x n illustrates the high-frequency part in the geomagnetic north component of the high-resolution map. Rather than considering the sparse codes for different channels independently, they can be jointly determined by reformulating Eq. (7) as Eq. (10).
where the cost function can be formulated as Eq. (11).
For simplicity, the regularization parameters λ and γ for all channels are set to be the same. The high-pass filters S n , S e and S v are also assumed to be the same for each channel. It should be noted that Eq. (11) is reduced to three independent sparse representation problems when τ = 0. In contrast to the single-channel sparse representation problem, the sparse codes in Eq. (10) are no longer independent. Therefore, Eq. (11) presents a more challenging optimization problem. Next, with the inspiration of [12], [25], a compliant solution for the new problem is proposed.
The following matrices and vectors are defined: where α c and y c constitute a concatenation of sparse codes and LR patches, respectively, in the three channels. P and P s are shifting matrices. The sparse code length of each channel is m, and the size of the HR patches is p. D l c and D h c are dictionaries that include the dictionaries of the three channels in their block diagonals. To simplify the formulation, we define D hs as Eq. (12).
and the cost function can be rewritten as Eq. (13).
and the optimization of α c is written as Eq. (15).
The problem in Eq. (15) is a convex sparsity constrained optimization, which can be solved with the algorithm in [26].

C. DICTIONARY INITIALIZATION
Before dictionary training, prior geomagnetic information is used for dictionary initialization. Since the collected geomagnetic field information is not equally spaced [7], this paper first models the geomagnetic field of the target area and then grids according to the modeling results represented by Eq. .
where B R is the residual geomagnetic field value, B 0 is the observed value for the geomagnetic field, and B C is the theoretical value calculated from the international geomagnetic reference field (IGRF). The residual magnetic field values of the target areas can be determined from Eq. (16). In a prior study, a model of the geomagnetic field is obtained by using spherical harmonic analysis (SHA) and rectangular harmonic analysis (RHA). In this paper, geomagnetic reference maps are determined from small areas where the effect of the earth's curvature is negligible. Therefore, RHA is employed to build a partial geomagnetic model based on B R . The actual procedure is shown below.
The research object of RHA is a rectangular region. In a space without a magnetic field source, the magnetic position satisfies the Laplace equation (Eq. (17)).
In Eq. (18), x, y, and z are the three-dimensional coordinates of the target points in the rectangular coordinate system, V (x, y, z) is the maximum truncation order, and L x , L y are the north-south and east-west length of the rectangular area. The center of the target area is selected as the origin of the rectangular harmonic coordinate system. The coordinates are restricted by: L y 2 and the three components of the magnetic field can be expressed by Eq. (19).
With the method above, a geomagnetic model can be constructed, from which we can derive the concrete gridded data. The initial redundant dictionary consists of 512 feature vectors. Half of these vectors are extracted from the international geomagnetic reference field, and the others are composed of the residual geomagnetic field modeled above.

D. DICTIONARY TRAINING
Redundant dictionary training is of considerable significance for achieving the desired results for sparse representation. By training the dictionaries of the LR and HR reference patches, same sparse representation can be achieved for LR and HR reference maps [27]. Sparse dictionaries are generally divided into two types: structure-based and learning-based. Structure-based dictionaries include Gaussian random matrices and circulant matrices. The most significant advantage of these dictionaries is that the calculation speed is very fast (with the complexity of o(n log n)) [28], but the ability to sparse represent the target signals is limited. Also, these types of dictionaries generally only prove useful on certain types of signals. A learningbased approach learns from a pre-constructed training dataset to build a redundant dictionary rather than from some specialized model.
In [23], the author applied the K-SVD method for image denoising and achieved good results. In this study, we used the K-SVD method to learn low-resolution and highresolution sparse dictionaries. The Frobenius norm is applied to measure the deviation between the original image and the sparsely decomposed image. The optimal problem is formulated in Eq. (21).
where d j stands for the j-th column of D and x T j stands for the j-th of X , keeping all the columns in the sparse matrix D except j 0 -th fixed. After extracting d j 0 individually, Eq. (21) is rewritten as Eq. (22).
The term in parentheses in Eq. (22) is defined as Eq. (23).
The minimization problem (Eq. (22)) is equal to the approximation of E j 0 , which can be solved via the singular value decomposition (SVD). However, this may generate a dense vector x T j , indicating that the non-zeros in X increase. The columns in E j 0 should be extracted -those columns where the corresponding entries in x T j 0 are non-zero -to keep the sparsity of X sufficient.
Therefore, a restriction matrix P j 0 is defined, which multiplies E j 0 from the right to preserve the relevant columns. Similarly, we define (x R j 0 ) T = x T j 0 P j 0 to preserve the non-zeros of x T j 0 only, and then Eq. (22) can be rewritten as Eq. (24).
An approximation via SVD can be applied to solve Eq. (24), by alternatively updating d j 0 and the corresponding coefficients in x T j 0 . The concrete solution is shown as follows: 1)Keeping d j 0 fixed and updating x R j 0 by solvinĝ 2)Keeping x R j 0 fixed and updating d j 0 once Eq. (24) is updated by solvingd The desired dictionary can be achieved after several iterations. Fig. 5 shows the K-SVD method in detail.

IV. EXPERIMENTS A. EXPERIMENT SETTINGS
Experiments of the proposed method are performed on a personal notebook ThinkPad T450, which has a configuration of an Intel Core i5-6200U CPU @ 2.30 GHz, 4 GB memory. MATLAB 2014a is used for simulation and experimental analysis. The database used in our experiments is a combination of two parts. One part is based on the geomagnetic anomaly data in North America issued by the United States Geological Survey (USGS) in 2002 [24]. The second part is the main magnetic field strength calculated from IGRF. We compare the proposed SR-GRM method with several well-known geomagnetic reference maps and single image super-resolution methods. These include bicubic interpolation, Kriging interpolation, and PSO-Kriging interpolation because they are all widely used in geomagnetic reference map construction. Another method we report results is Neighbor Embedding (NE). First, the values of the geomagnetic reference maps are normalized to the [0,255] range linearly, according VOLUME 8, 2020 to Eq. (27), where t i denotes the pixel in the geomagnetic reference maps, Max denotes the maximum value, and Min denotes the minimum value in all maps.
In our experiments, the input geomagnetic reference maps are magnified by a factor of 2, 3, or 4. For the LR maps, we use 5×5 LR patches with an overlap of 4 pixels between adjacent patches and redundant dictionaries based on the algorithm in [29].  The features are extracted from the bicubic interpolated version of the whole map with the target magnification factor rather than from the 5 × 5 LR patch directly. Then, extracted features are applied to generate the sparse coefficients according to [15], [29], [30]. Finally, the HR patches are reconstructed from the same sparse coefficients and the HR dictionary. Dictionaries are trained in over 50000 patch pairs in the geomagnetic database. The number of features included in each trained dictionary is 512 and λ is set via cross-validation of 0.2. Five hundred reference maps are extracted randomly as the validation set.

B. EXPERIMENTAL RESULTS
In order to evaluate the obtained super-resolution reference maps quantitatively, five frequently-used image quality matrices are applied, which are the Peak Signal to Noise Ratio (PSNR), Structural Similarity Index (SSIM), Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Standard Deviation (SD). The total geomagnetic field strength is calculated for quality evaluation using Eq. (20).
In Fig. 6, we compare the reconstruction result of our algorithm with several other methods, including neighbor embedding, bicubic, and PSO-Kriging. The magnification factor is 4X, and the size of the geomagnetic reference map is 128 × 128. PSO-Kriging interpolation and our method give the best results. Our proposed method better highlights the details of the reference map.
To verify the superiority of multi-channel sparse representation method, single-channel method is used for comparison. The scaling factor is 3. The experimental results are shown in Table 1. Multi-channel method outperforms single-channel method in both PSNR, SSIM and RMSE indicators. Tables 2-4 show the super-resolution reconstruction results of maps in the validation set. PSNR, SSIM, and RMSE are compared, and our method outperforms the other competing methods. The advantage of our method increases with the magnification factor. Bicubic interpolation is the most commonly used method for its low computational complexity.  However, its performance in super-resolution reconstruction decreases rapidly with large magnification factors. Neighbor embedding generates sharp edges in places, but the textures are blurred among the pixels of LR maps. PSO-Kriging somewhat adapts to the geomagnetic reference map reconstruction but generates undesired smoothing that is not present in the original HR maps. Learning-based methods such as Generative Adversarial Network (GAN) and Convolutional Neural Network (CNN) are also added for comparison. The training set size is the same as that used in the proposed method. Results in Tables 2-4 indicate that the performance of learning-based methods are unsatisfying while the training set is not sufficiently big. So far, we have used a fixed training set of size 50000 for all the experiments. In the following experiment, we verify the importance of prior information for dictionary initialization, which has been discussed in part B of section III. We again trained dictionaries of size 512 atoms and samples of 1600, 3200, 6400, 12800, 25600, 51200 map patches as training sets. The results are also evaluated in terms of PSNR, SSIM, and RMSE. Figures 7-9 show the variation of dictionaries initialized with prior information and random entries versus training set size. When prior information is attached, the dictionary can reach the same reconstruction precision with a relatively small training set size.   An often-made assumption is that the input reference maps are free of noise, which is violated in real geomagnetic data collections. The artifacts introduced in the denoising process may remain and are even magnified after the super-resolution   procedure. Different levels of Gaussian noise are attached to the LR reference maps to test the robustness of our method to noise. Similar to [29], the range of the standard deviation of the noise is chosen from 3 to 12 with a magnification factor of 3. Table 5-7 shows the PSNR, SSIM, and RMSE results of reconstruction maps with different levels of noise. The proposed method outperforms the competition in all cases.

V. CONCLUSION
In this paper, we proposed a novel method for geomagnetic reference map reconstruction that makes a considerable improvement over existing reconstruction methods. Besides producing high precision reconstruction results, our method can significantly reduce the impact of noise, which contributes to actual geomagnetic reference map construction. Furthermore, a dictionary initialization method based on prior information of RHA is proposed to train a more qualitative sparse dictionary with a relatively small training set size. In future work, we will investigate the coherence of the three independent components of the geomagnetic reference field and further improve the precision of the geomagnetic reference map.