地磁匹配导航技术主要使用实时采集的地磁场数据与该区域地磁基准图进行匹配[3-4]，从而实现定位。因此地磁基准图的构建精确度对定位精度的好坏起着决定性作用。目前构建地磁基准图的方法主要有2种：一是根据现有的地磁场物理模型进行构建; 二是根据实测地磁场数据，构建网格化地磁基准图[5]。现有国际地磁参考磁场(IGRF)和世界磁场模型(WMM)主要描述的是地球主磁场模型[6]。地磁场模型依赖卫星监测数据，地球内部异常场随海拔高度升高而衰减，因此地磁场模型中几乎不含异常场信息，用于刻画局部地磁场时误差很大，通过实测数据获得的地磁基准图具有更大优势。目前，地磁基准图的构建主要集中于插值法的研究，常用的方法有：三次样条插值、Kriging插值[7]、基于粒子群优化(PSO)算法的Kriging(PSO-Kriging)插值[8]等。
现有插值方法存在的共同问题是在实测数据量较小时误差较大，而压缩感知理论能够在采集少量数据的同时几乎完全重构出原始信号，可以弥补在低采样率下插值方法的缺陷。因此本文以现有的地磁场实测数据为基础，利用信号处理中的压缩感知理论进行地磁基准图的构建，研究适用于地磁信息特征的压缩感知算法，并验证其在实测数据量较小时相比于传统插值法的优越性。
1 压缩感知理论 压缩感知是2006年由Candes[9]和Donoho[10-11]等提出的具有革命性的信号处理理论。压缩感知理论指出，当信号在某个变换域上是稀疏的，就可以投影为低维观测向量，以远低于奈奎斯特理论的采样频率进行信号采集，同时能够高精度甚至完全重构出原始信号[12]。其核心思想是在信号发送端进行压缩采样，并在信号接收端进行高精度重建的过程，目的是在保证信号精度的同时大幅度提升信号采集的效率并降低采集成本[13]。压缩感知理论主要分为稀疏表示和信号重建2个部分。设一个信号x(x∈RN×1)，在基底Ψ下表示为信号α，即x=Ψα。若α中只有K个元素非零或其他元素趋近于零，则称该信号在基底Ψ下是K稀疏的，Ψ称为稀疏字典。
若能找到合适的观测矩阵Φ∈RM×N(K <
 | (1) |
式中：A为传感矩阵。即已知观测信号y，求解出稀疏信号α，从而得到初始信号x。
观测信号的维度M大于α的稀疏度K，因此可证明存在唯一解，即在已知压缩矩阵Ψ和观测矩阵Φ时，根据低维信号y必然可解出原始高维信号x。因此实际上该问题可转化为在一定稀疏度条件下求解最小L0范数解的问题，即
 | (2) |
式中：b为原始信号；A=ΦΨ。显然Ψ和Φ共同决定了变换后信号的稀疏性，Ψ与Φ的差异越大，得到的解α就越稀疏。因此定义Ψ与Φ的互相关性为
 | (3) |
其中：
要想将原本n维的信号降到m维同时保留绝大部分重要信息，根据压缩感知理论，传感矩阵A需要满足约束等距性(Restricted Isometry Property，RIP)[14]条件，即传感矩阵A中的RIP参数δK为满足下式的最小值：
 | (4) |
式中：若δK < 1，则称传感矩阵A满足K阶RIP条件。
由于求L0范数解是一个NP(Non-deterministic Polynomial)难问题，在维度较高时难以求解；在满足RIP条件时，L0范数与L1范数等价，从而可转化为一个求最大L1范数的凸优化问题进行求解，即
 | (5) |
观测矩阵Φ的维度直接决定了压缩感知的效果，过大就失去压缩的意义，过小会导致重要信息的丢失，满足：
 | (6) |
式中：C为常数; m为压缩观测后信号的维度；μ(A)为传感矩阵互相关系数；r为信号的稀疏度。由于C的值难以确定以及μ(A)计算困难，一般使用经验公式(7)计算：
 | (7) |
2 稀疏字典 现有常用的稀疏字典有离散傅里叶基、离散余弦基、离散小波变换基、Curvelet基、冗余字典等[15]。
2.1 离散傅里叶变换 离散傅里叶变换将二维信号利用一系列周期函数进行表示，实际上是从时域转化到频域进行分析，离散傅里叶基可表示为
 | (8) |
式中：n由原始信号的长度决定。若以离散傅里叶基作为稀疏基，单位阵In×n作为观测矩阵，则根据式(3)可得
2.2 离散余弦变换 离散余弦变换是原信号在进行偶延拓后展开成的仅含余弦项的傅里叶级数，离散余弦变换矩阵为
 | (9) |
若以单位阵In×n作为观测矩阵，根据式(3)可得
2.3 离散小波变换 离散小波变换是对于连续小波的尺度和平移进行离散化，在图像处理方面有显著的应用效果。本文采用Haar小波变换作为小波基。二维小波变换过程如图 1所示。
 |
| 图 1 二维小波变换过程 Fig. 1 Two-dimensional wavelet transform process |
| 图选项 |
3 观测矩阵 3.1 高斯随机矩阵 高斯随机矩阵是在压缩感知中应用最广泛的一种观测矩阵，由于其具有很强的随机性，只需观测维度M满足：
 | (10) |
式中:c为常数即可高概率地满足RIP条件。因此在压缩感知算法仿真中一般使用高斯随机矩阵作为观测矩阵以达到较好的实验效果。但由于构造矩阵时数据的连续性和随机性，硬件上实现困难，因此难以得到实际应用。
3.2 随机伯努利矩阵 随机伯努利矩阵的性质类似于高斯随机矩阵，由于其中每一个元素都服从伯努利分布，M≥cKlg(N/K)时以接近1的概率满足RIP条件。相比于高斯随机矩阵，由于其元素都为±1，更有利于硬件的实现。
3.3 稀疏随机矩阵 稀疏矩阵构造方法为：在一个M×N的全零阵基础上，将每一行随机d个位置置1，其中d∈{4, 8, 10, 16}，在保证一定的重建精度的条件下易于进行构造和实际实现。
3.4 循环矩阵 循环矩阵的构造方法为：首先构造随机向量β，即β=(β1, β2, …, βN)∈RN，β中的元素取±1，而后经过M次移位操作构造出完整的矩阵，同时矩阵中每一个元素服从伯努利分布。由于可通过移位寄存器实现矩阵构造，具有较好的应用前景。
4 仿真结果与性能分析 为研究适用于地磁基准图构建的压缩感知方法，本文比较了多种稀疏基和观测矩阵在地磁基准图构建中的效果。稀疏基选择离散小波基和离散余弦基，观测矩阵选择高斯随机矩阵、随机伯努利矩阵、稀疏随机矩阵、循环矩阵以及单位矩阵。本文以美国地质调查局(United States Geological Survey，USGS)于2002年发布的北美地区地磁场异常数据[16]为基础，结合IGRF计算得到的地磁主磁场强度，从而得到地磁总磁场实验数据。实验样本为10张128×128网格具有不同起伏程度的地磁基准图，选取区域大小都为约50 km×50 km，以保证基本涵盖全部类型的地磁特征；使用观测矩阵维度、图像处理中的峰值信噪比(Peak Signal to Noise Ratio，PSNR)以及绝对误差(MAE)和均方根误差(RMSE)作为评价标准。由于在信号处理端无需考虑实时性问题，信号重构算法采用精度较高的压缩采样匹配追踪(CoSaMP)算法，通过不断更新迭代原子提高精度。PSNR值计算公式为
 | (11) |
 | (12) |
式中：I为原图像；G为处理图像。
4.1 稀疏矩阵分析 以高斯随机矩阵作为观测矩阵，检验离散小波变换和离散余弦变换在不同维度观测矩阵下的重建效果。图 2为离散小波变换和离散余弦变换重建效果随观测矩阵维度变化曲线，PSNR值为10张样本图分别重构后的平均值。可以看出，当观测矩阵维度较大时离散小波变换重建效果较高，但随着观测矩阵维度的降低，其重建性能迅速降低，M < 50时几乎无法重建信号；而对于离散余弦变换，M对重建效果影响较小，即使降到70，PSNR值依然可以保持在80 dB以上。
 |
| 图 2 离散小波变换和离散余弦变换重建效果随观测矩阵维度的变化 Fig. 2 Variation of reconstruction effect of discrete wavelet transform and discrete cosine transform with dimension of measurement matrix |
| 图选项 |
由于取的信号稀疏度K=M/4，当M减小时信号稀疏度也随之减小。经过离散小波变换后的信号难以保证足够的稀疏度，因此当观测矩阵维度降低时会丢失大量的信息，造成重建信号的严重失真；而离散余弦变换可以保证足够的压缩率，能够在保证一定重建精度的条件下降低观测信号的维度。
4.2 观测矩阵分析 为寻找适合地磁基准图的观测矩阵，以离散余弦变换作为稀疏矩阵，分别以高斯随机矩阵、随机伯努利矩阵、稀疏随机矩阵、循环矩阵、单位矩阵作为观测矩阵检验重建效果，实验效果如图 3所示。
 |
| 图 3 不同观测矩阵重建效果随观测矩阵维度的变化 Fig. 3 Variation of reconstruction effect with dimension of matrix for different measurement matrices |
| 图选项 |
从图 3可以看出，当观测矩阵维度较大时，循环矩阵、随机伯努利矩阵和高斯随机矩阵都具有较好的信号重构效果; 但随着维度的降低，高斯随机矩阵、随机伯努利矩阵、稀疏随机矩阵、循环矩阵的重构性能都迅速降低; 在M < 55时，单位矩阵由于性能基本不受观测矩阵维度影响逐渐体现出其明显优势。因此在保证一定重构精度条件下使用单位矩阵作为观测矩阵能够减小观测信号的维度。
为检验重建算法的稳定性，将采样率设为30%，对每种观测矩阵进行100次重建，重建效果如图 4所示。
 |
| 图 4 不同观测矩阵重建稳定性 Fig. 4 Reconstruction stability of different measurement matrices |
| 图选项 |
以PSNR值的标准差(SD)作为判定重建算法稳定性的标准，计算结果如表 1所示，在30%采样率下，利用单位矩阵作为观测矩阵，在重建精度和稳定性都优于其他观测矩阵。
表 1 不同观测矩阵重建100次PSNR标准差 Table 1 PSNR standard deviation of 100 times reconstruction for different measurement matrices ?
| 观测矩阵 | PSNR标准差 |
| 高斯随机矩阵 | 4.9845 |
| 随机伯努利矩阵 | 4.8488 |
| 稀疏随机矩阵 | 5.3112 |
| 循环矩阵 | 4.5170 |
| 单位矩阵 | 0 |
表选项
4.3 与其他重建算法比较 为检验CoSaMP算法的效果，分别使用匹配追踪(MP)算法、正交匹配追踪(OMP)算法、分段正交匹配追踪(StOMP)算法对30%采样率下的地磁基准图进行重建，以单位矩阵作为观测矩阵，以PSNR值作为评价标准，实验结果如表 2所示。
表 2 不同重建算法结果 Table 2 Results of different reconstruction algorithms ?
| 重建算法 | PSNR值 |
| CoSaMP | 91.45 |
| MP | 88.33 |
| OMP | 90.12 |
| StOMP | 90.66 |
表选项
4.4 压缩感知与插值法比较 为验证压缩感知在地磁基准图构建中的实际效果，将离散余弦变换及单位矩阵构成的传感矩阵和三次样条插值、Kriging插值、PSO-Kriging插值分别在25%、18.75%、12.5%、6.25%采样率下的基准图构建效果进行对比。
Kriging插值所采用的变异函数模型为指数模型，其具体形式为
 | (13) |
式中：T0为基台值；T为块金值；h为到插值点的距离；3a为变程。PSO-Kriging插值变异函数同样采用指数模型，PSO算法粒子总个数设置40个，迭代次数设置100次。
具体采样方法为：设128×128的原始地磁基准图为P=pij(i, j=1, 2, …, 128)，在25%采样率(M/N)时，压缩感知传感矩阵维度选择32×128，即M=32在6.25%采样率时，传感矩阵维度选择8×128，插值所用基准图由随机下采样获得，即保证所有方法采样数据量相同。以PSNR值作为评价指标，结果如表 3所示。当采样率为25%时，压缩感知与插值法的PSNR值基本相近，压缩感知精度略高于插值法；当采样率降到6.25%时，压缩感知依然可以保持70 dB以上的PSNR值，明显优于同采样率下的插值法。
表 3 不同采样率下各方案构建基准图的PSNR值 Table 3 PSNR values of reference map reconstructed by different schemes at different sampling rates
| 采样率/% | PSNR/dB | |||
| 压缩感知 | 三次样条插值 | Kriging插值 | PSO-Kriging插值 | |
| 25 | 88.2135 | 84.5327 | 86.4336 | 86.8971 |
| 18.75 | 83.4127 | 79.1009 | 81.3121 | 81.6743 |
| 12.5 | 78.3155 | 70.4422 | 72.4561 | 74.7963 |
| 6.25 | 74.6714 | 60.5773 | 64.3352 | 66.9712 |
表选项
图 5为用压缩感知和Kriging插值重建后的基准图效果，可以看出Kriging插值后失去了很多细节信息，经过压缩感知后图像与原始图相似度更高。
 |
| 图 5 Kriging插值法和压缩感知重构效果对比 Fig. 5 Comparison of Kriging interpolation method and compressed sensing reconstruction effect |
| 图选项 |
为充分验证实验结果的有效性，利用文献[5]中的绝对误差、均方根误差和标准差3种性能指标进行分析。
绝对误差为
 | (14) |
均方根误差为
 | (15) |
标准差为
 | (16) |
式中：ξi为真值与估计值之间的误差。
为验证各种方法构建效果随采样率变化情况，本文选取了25%、18.75%、12.5%、6.25%这4种采样率进行实验，计算结果分别如表 4~表 7所示。
表 4 采样率下各性能指标数据 Table 4 Performance index data at 25% sampling rate ?
| 构建方法 | 绝对误差 | 均方根误差 | 标准差 |
| 压缩感知 | 0.5978 | 0.6119 | 0.6743 |
| 三次样条插值 | 2.4533 | 3.2121 | 3.2156 |
| Kriging插值 | 2.2002 | 3.0117 | 3.1189 |
| PSO-Kriging插值 | 2.0524 | 2.3768 | 2.5143 |
表选项
表 5 18.75%采样率下各性能指标数据 Table 5 Performance index data at 18.75% sampling rate ?
| 构建方法 | 绝对误差 | 均方根误差 | 标准差 |
| 压缩感知 | 3.1033 | 3.4561 | 3.4788 |
| 三次样条插值 | 6.2455 | 7.5371 | 7.6133 |
| Kriging插值 | 5.6413 | 6.4745 | 6.2531 |
| PSO-Kriging插值 | 5.0177 | 5.7874 | 5.8125 |
表选项
表 6 12.5%采样率下各性能指标数据 Table 6 Performance index data at 12.5% sampling rate ?
| 构建方法 | 绝对误差 | 均方根误差 | 标准差 |
| 压缩感知 | 5.4177 | 6.0133 | 6.3145 |
| 三次样条插值 | 18.7133 | 24.9103 | 25.3174 |
| Kriging插值 | 11.9773 | 15.6245 | 15.6711 |
| PSO-Kriging插值 | 11.3144 | 14.5178 | 14.6477 |
表选项
表 7 6.25%采样率下各性能指标数据 Table 7 Performance index data at 6.25% sampling rate ?
| 构建方法 | 绝对误差 | 均方根误差 | 标准差 |
| 压缩感知 | 10.2563 | 11.3313 | 11.8847 |
| 三次样条插值 | 37.3318 | 41.6547 | 41.7734 |
| Kriging插值 | 26.8891 | 30.3671 | 31.2155 |
| PSO-Kriging插值 | 25.4731 | 28.5717 | 29.0973 |
表选项
从表 4~表 7可以看出，在各种采样率下利用压缩感知理论重构地磁基准图性能都优于插值法，且采样率越低优势越明显。
5 结论 1) 验证了压缩感知理论在地磁基准图构建中的可行性以及稳定性。
2) 确定了以离散余弦变换作为稀疏基、单位矩阵作为观测矩阵，以CoSaMP算法作为重构算法的基于压缩感知的高精度地磁基准图构建方法。
3) 仿真结果表明，在6.25%采样率下通过压缩感知理论重构的地磁基准图PSNR值、平均误差、均方根误差指标都明显优于插值方法。
4) 根据压缩感知理论，要想完成对信号的压缩感知，首先要进行压缩采样，通过一个随机矩阵控制采样点选择是否采集信号。这就需要重新设计一种地磁信息采集方案。
5) 如何将压缩感知理论实际应用到地磁传感器设计中并设计低成本、高重构精度同时利于硬件实现的观测矩阵将是今后的研究方向。
参考文献
| [1] | GOLDENBERG F. Geomagnetic navigation beyond the magnetic compass[C]//Position, Location, & Navigation Symposium.Piscataway, NJ: IEEE Press, 2006: 684-694. |
| [2] | FOX S M, PAL P K, PSIAKI M L.Magnetometer-based autonomous satellite navigation(MAGNAV)[C]//Proceedings of the Annual Rocky Mountain Guidance and Control Conference.Washington, D.C.: American Astronautical Society, 1990: 369-382. |
| [3] | 蔡洪, 郭才发, 胡正东. 惯性/地磁组合导航算法[J]. 中国惯性技术学报, 2009, 17(3): 333-337. CAI H, GUO C F, HU Z D. Algorithms for inertial/geomagnetic integrated navigation[J]. Journal of Chinese Inertial Technology, 2009, 17(3): 333-337. (in Chinese) |
| [4] | CHEON Y J. Fast convergence of orbit determination using geomagnetic field measurement in target pointing satellite[J]. Aerospace Science and Technology, 2013, 30(1): 315-322. DOI:10.1016/j.ast.2013.08.016 |
| [5] | 王哲, 王仕成, 张金生, 等. 一种地磁匹配制导基准图构建方法及其有效性评价[J]. 系统工程与电子技术, 2008, 30(11): 2207-2211. WANG Z, WANG S C, ZHANG J S, et al. Method for preparation of reference map in geomagnetism matching guidance and its validity evaluation[J]. Systems Engineering and Electronics, 2008, 30(11): 2207-2211. DOI:10.3321/j.issn:1001-506X.2008.11.041 (in Chinese) |
| [6] | 踪华, 刘嬿, 杨业. 地磁导航技术研究现状综述[J]. 航天控制, 2018, 36(3): 94-99. ZONG H, LIU Y, YANG Y. Overview of the research status about geomagnetic navigation technology[J]. Aerospace Control, 2018, 36(3): 94-99. (in Chinese) |
| [7] | 岳建平, 甄宗坤. 基于粒子群算法的Kriging插值在区域地面沉降中的应用[J]. 测绘通报, 2012(3): 59-62. YUE J P, ZHEN Z K. Application of particle swarm optimization based Kriging interpolation method in regional land subsidence[J]. Bulletin of Surveying and Mapping, 2012(3): 59-62. (in Chinese) |
| [8] | 李晨霖, 王仕成, 张金生, 等. 基于改进的Kriging插值方法构建地磁基准图[J]. 计算机仿真, 2018, 35(12): 278-282. LI C L, WANG S C, ZHANG J S, et al. Construction of geomagnetic datum map based on improved Kriging interpolation method[J]. Computer Simulation, 2018, 35(12): 278-282. (in Chinese) |
| [9] | CANDES E J, ROMBERG J, TAO T. Robust uncertainty principles:Exact signal reconstruction from highly incomplete frequency information[J]. IEEE Transactions on Information Theory, 2006, 52(2): 489-509. DOI:10.1109/TIT.2005.862083 |
| [10] | DONOHO D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306. DOI:10.1109/TIT.2006.871582 |
| [11] | DONOHO D L, TSAIGY. Extensions of compressed sensing[J]. Signal Processing, 2006, 86(3): 533-548. DOI:10.1016/j.sigpro.2005.05.028 |
| [12] | FOUCARTS, RAUHUT H.A mathematical introduction to compressive sensing[M].[S.l.]: Birkhauser, 2013: 31-37. |
| [13] | THEODORIDIS S, KOPSINIS Y, SLAVAKIS K. Sparsity-aware learning and compressed sensing:Anoverview[J]. Mathematics, 2012, 1: 1271-1377. |
| [14] | MOHADES M M, MOHADES A, TADAION A. A reed-solomon code based measurement matrix with small coherence[J]. IEEE Signal Processing Letters, 2014, 21(7): 839-843. DOI:10.1109/LSP.2014.2314281 |
| [15] | 吴赟.压缩感知测量矩阵的研究[D].西安: 西安电子科技大学, 2012: 32-46. WU Y.Research on measurement matrix for compressive sensing[D].Xi'an: Xidian University, 2012: 32-46(in Chinese). |
| [16] | U.S.Geological Survey.Digital data grids for the magnetic anomaly map of North America: open-file report 02-414[R].Denver: Geological Survey, 2002. |
