马啸宇，张金生，*，郝亮亮，李婷，王少博，李琳

（1.火箭军工程大学 精确制导与仿真实验室，西安710025； 2.火箭军96863部队，洛阳471000）

在近年来的导航技术发展过程中，卫星导航和惯性导航技术一直是飞行器导航技术的主流发展方向。但卫星导航的抗干扰能力较差，惯性导航的误差随时间不断累积。地磁场随时间变化缓慢且不受气候、地域影响，且由于其无需主动发射探测信号、抗干扰能力强，地磁导航技术作为一种辅助导航方式逐渐受到关注［1-2］。

地磁匹配导航技术主要使用实时采集的地磁场数据与该区域地磁基准图进行匹配［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＜M≪N）对信号x进行采样，则可得到M 维观测信号为

式中：A为传感矩阵。即已知观测信号y，求解出稀疏信号α，从而得到初始信号x。

观测信号的维度M大于α的稀疏度K，因此可证明存在唯一解，即在已知压缩矩阵Ψ 和观测矩阵Φ时，根据低维信号y必然可解出原始高维信号x。因此实际上该问题可转化为在一定稀疏度条件下求解最小L0范数解的问题，即

式中：b为原始信号；A＝ΦΨ。显然Ψ 和Φ 共同决定了变换后信号的稀疏性，Ψ 与Φ 的差异越大，得到的解α就越稀疏。因此定义Ψ 与Φ 的互相关性为

2 稀疏字典

现有常用的稀疏字典有离散傅里叶基、离散余弦基、离散小波变换基、Curvelet基、冗余字典等［15］。

2.1 离散傅里叶变换

离散傅里叶变换将二维信号利用一系列周期函数进行表示，实际上是从时域转化到频域进行分析，离散傅里叶基可表示为

2.2 离散余弦变换

2.3 离散小波变换

离散小波变换是对于连续小波的尺度和平移进行离散化，在图像处理方面有显著的应用效果。本文采用Haar小波变换作为小波基。二维小波变换过程如图1所示。

图1 二维小波变换过程Fig.1 Two-dimensional wavelet transform process

3 观测矩阵

3.1 高斯随机矩阵

高斯随机矩阵是在压缩感知中应用最广泛的一种观测矩阵，由于其具有很强的随机性，只需观测维度M满足：

式中：c为常数即可高概率地满足RIP条件。因此在压缩感知算法仿真中一般使用高斯随机矩阵作为观测矩阵以达到较好的实验效果。但由于构造矩阵时数据的连续性和随机性，硬件上实现困难，因此难以得到实际应用。

3.2 随机伯努利矩阵

随机伯努利矩阵的性质类似于高斯随机矩阵，由于其中每一个元素都服从伯努利分布，M≥cK lg（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）作为评价标准。由于在信号处理端无需考虑实时性问题，信号重构算法采用精度较高的压缩采样匹配追踪（Co-SaMP）算法，通过不断更新迭代原子提高精度。PSNR值计算公式为

4.1 稀疏矩阵分析

以高斯随机矩阵作为观测矩阵，检验离散小波变换和离散余弦变换在不同维度观测矩阵下的重建效果。图2为离散小波变换和离散余弦变换重建效果随观测矩阵维度变化曲线，PSNR值为10张样本图分别重构后的平均值。可以看出，当观测矩阵维度较大时离散小波变换重建效果较高，但随着观测矩阵维度的降低，其重建性能迅速降低，M＜50时几乎无法重建信号；而对于离散余弦变换，M对重建效果影响较小，即使降到70，PSNR值依然可以保持在80 dB以上。

由于取的信号稀疏度K＝M／4，当M 减小时信号稀疏度也随之减小。经过离散小波变换后的信号难以保证足够的稀疏度，因此当观测矩阵维度降低时会丢失大量的信息，造成重建信号的严重失真；而离散余弦变换可以保证足够的压缩率，能够在保证一定重建精度的条件下降低观测信号的维度。

图2 离散小波变换和离散余弦变换重建效果随观测矩阵维度的变化Fig.2 Variation of reconstruction effect of discrete wavelet transform and discrete cosine transform with dimension of measurementmatrix

4.2 观测矩阵分析

为寻找适合地磁基准图的观测矩阵，以离散余弦变换作为稀疏矩阵，分别以高斯随机矩阵、随机伯努利矩阵、稀疏随机矩阵、循环矩阵、单位矩阵作为观测矩阵检验重建效果，实验效果如图3所示。

图3 不同观测矩阵重建效果随观测矩阵维度的变化Fig.3 Variation of reconstruction effect with dimension of matrix for differentmeasurementmatrices

从图3可以看出，当观测矩阵维度较大时，循环矩阵、随机伯努利矩阵和高斯随机矩阵都具有较好的信号重构效果；但随着维度的降低，高斯随机矩阵、随机伯努利矩阵、稀疏随机矩阵、循环矩阵的重构性能都迅速降低；在M＜55时，单位矩阵由于性能基本不受观测矩阵维度影响逐渐体现出其明显优势。因此在保证一定重构精度条件下使用单位矩阵作为观测矩阵能够减小观测信号的维度。

为检验重建算法的稳定性，将采样率设为30%，对每种观测矩阵进行100次重建，重建效果如图4所示。

以PSNR值的标准差（SD）作为判定重建算法稳定性的标准，计算结果如表1所示，在30%采样率下，利用单位矩阵作为观测矩阵，在重建精度和稳定性都优于其他观测矩阵。

图4 不同观测矩阵重建稳定性Fig.4 Reconstruction stability of different measurementmatrices

表1 不同观测矩阵重建100次PSNR标准差Table 1 PSNR standard deviation of 100 tim es reconstruction for d ifferentm easurem ent m atrices dB

4.3 与其他重建算法比较

为检验CoSaMP算法的效果，分别使用匹配追踪（MP）算法、正交匹配追踪（OMP）算法、分段正交匹配追踪（StOMP）算法对30%采样率下的地磁基准图进行重建，以单位矩阵作为观测矩阵，以PSNR值作为评价标准，实验结果如表2所示。

4.4 压缩感知与插值法比较

为验证压缩感知在地磁基准图构建中的实际效果，将离散余弦变换及单位矩阵构成的传感矩阵和三次样条插值、Kriging插值、PSO-Kriging插值分别在25%、18.75%、12.5%、6.25%采样率下的基准图构建效果进行对比。

Kriging插值所采用的变异函数模型为指数模型，其具体形式为

式中：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值，明显优于同采样率下的插值法。

图5为用压缩感知和Kriging插值重建后的基准图效果，可以看出Kriging插值后失去了很多细节信息，经过压缩感知后图像与原始图相似度更高。

为充分验证实验结果的有效性，利用文献［5］中的绝对误差、均方根误差和标准差3种性能指标进行分析。

绝对误差为

表2 不同重建算法结果Tab le 2 Results of differen t reconstruction algorithm sdB

表3 不同采样率下各方案构建基准图的PSNR值Tab le 3 PSNR values of reference m ap reconstructed by different schem es at differen t sam p ling rates

图5 Kriging插值法和压缩感知重构效果对比Fig.5 Comparison of Kriging interpolation method and compressed sensing reconstruction effect

均方根误差为

式中：ξi为真值与估计值之间的误差。

为验证各种方法构建效果随采样率变化情况，本文选取了25%、18.75%、12.5%、6.25%这4种采样率进行实验，计算结果分别如表4～表7所示。

从表4～表7可以看出，在各种采样率下利用压缩感知理论重构地磁基准图性能都优于插值法，且采样率越低优势越明显。

表4 25%采样率下各性能指标数据Table 4 Perform ance index data at 25% sam p ling rate nT

表5 18.75%采样率下各性能指标数据Table 5 Perform ance index data at 18.75%sam p ling rate nT

表6 12.5%采样率下各性能指标数据Table 6 Perform ance index data at 12.5%sam p ling rate nT

表7 6.25%采样率下各性能指标数据Tab le 7 Perform ance index data at 6.25%sam p ling rate nT

5 结 论

1）验证了压缩感知理论在地磁基准图构建中的可行性以及稳定性。

2）确定了以离散余弦变换作为稀疏基、单位矩阵作为观测矩阵，以CoSaMP算法作为重构算法的基于压缩感知的高精度地磁基准图构建方法。

3）仿真结果表明，在6.25%采样率下通过压缩感知理论重构的地磁基准图PSNR值、平均误差、均方根误差指标都明显优于插值方法。

4）根据压缩感知理论，要想完成对信号的压缩感知，首先要进行压缩采样，通过一个随机矩阵控制采样点选择是否采集信号。这就需要重新设计一种地磁信息采集方案。

5）如何将压缩感知理论实际应用到地磁传感器设计中并设计低成本、高重构精度同时利于硬件实现的观测矩阵将是今后的研究方向。