1 引 言

随着航空地磁测量精度的提高以及高精度插值算法的应用，地磁数据的网格分辨率达到百米量级左右，大大提高了数据的空间分辨率，其丰富的纹理信息，更利于应用地磁及其梯度数据进行辅助导航^{[1, 2, 3, 4, 5]}。但由于数据量的剧增，对导航数据库构建提出了较高的要求，研究大容量地磁数据的压缩方法，降低高分辨率数据库的存储容量具有重要意义。

传统数据压缩方法是以牺牲数据的高频细节信息为代价，通常运用BDCT(blocking discrete cosine transform)编解码技术对数据进行分块变换处理^[6]。较大分块数据相关性较高，可以实现较大的压缩比，但大尺寸分块数据压缩时又会丢失数据的平稳性，引入振铃效应(ringing effect)噪声和块效应(block effect)边界误差^{[7, 8, 9]}。由此重构的数据高频轮廓特征信息损失严重或引入重构噪声，增大了地磁辅助导航系统匹配失败或是误匹配的概率^{[10, 11]}，因此不适用于大容量的导航数据压缩。

压缩感知(compressed sensing,CS)作为一种新兴的稀疏信号精确重构理论已成为核磁共振和合成孔径雷达成像等图像处理领域的重要方法^{[12, 13]}。由于地磁数据能量多集中在低频分量上，而压缩感知对稀疏信号的随机采样矩阵通常为局部傅里叶变换系数星形采样，这种采样矩阵破坏了低频分量系数之间的相关性，影响了数据重构精度^{[14, 15]}。

为解决大容量地磁数据压缩的同时，又能较好地保存细节轮廓信息的问题，本文提出一种基于压缩感知理论的数据分频压缩及联合重构方法。

压缩感知是一种新兴理论，它根据信号本身或在变换域中具有的稀疏特性，利用远少于Nyquist采样定理所需的采样数据就可以精确地重构原始数据，因此，该理论在各个领域中展开广泛研究及应用^[16]。

假设离散信号x∈R^N在小波变换下的稀疏度为K，则原始信号可用一组小波变换的稀疏基Ψ={ψ_i}_i=1^N线性表示为

点击浏览公式

根据压缩感知理论，当随机测量矩阵Φ与稀疏变换矩阵Ψ满足约束等距性(restricted isometry property,RIP)条件，即可利用信号x线性测量获得的观测值y=Φx=Au，通过解决以下l₁范数最小值问题来实现原始数据的精确重构

点击浏览公式

点击浏览公式

l₁范数定义为||w||₁=∑_i|w_i|，最小化||u||₁保证求解的稀疏性，约束||Au－y||₂²≤σ²使求解与最优解具有一致性，式(3)是在稀疏变换域Ψ中寻找最优稀疏解。

对于时域中绝大多数图像数据，其离散梯度具有稀疏性，全变差(total variation,TV)正则化模型能够比较合理地保持图像的轮廓结构和边缘特征^[17]。因此，全变差方法常作为稀疏变换，替代l₁范数应用在压缩感知方法中^[12]。文献^{[18, 19]}结合TV正则化和l₁范数的优点，给出同时在稀疏变换域以及稀疏梯度时域中寻找最优解的压缩感知模型

点击浏览公式

点击浏览公式

由于TV正则化和l₁范数都具有非平滑、不可微分的特性，传统优化算法如共轭梯度法和匹配追踪算法在解决以上优化问题时，需要耗费大量的计算时间^{[20, 21]}。近期一些快速算法如TwIST、FISTA、TVAL3和RecPF等的出现，使得在求解压缩感知问题时可以节省大量的时间^{[22, 23, 24, 25]}。

对于导航用地磁数据，压缩时需要保留精度较高的边缘、轮廓和纹理等局部高频特征，基于压缩感知理论的数据分频压缩及联合重构方法主要包含两部分内容：频域余弦变换的分频压缩和基于压缩感知的联合重构。

分频压缩方法的主要过程为：

(1) 以较大范围的高精度地磁数据作为整体进行二维离散余弦变换，以将数据进行多尺度频率分解。

(2) 设置变换域系数阈值α，将幅值大于阈值的低频分量系数进行无损编码。

(3) 高频分量系数包含边缘和纹理成分，利用放射线采样矩阵R对高频分量系数进行稀疏采样编码。

数据重构方法的主要过程为：

(1) 对分频压缩存储的低频系数和高频系数按所处位置进行重组，作为压缩感知重构中的线性观测值y=Φx。其中，测量矩阵Φ为局部余弦变换系数矩阵Φ=C_Ω=Ω′·C[·]，C[·]代表二维离散余弦变换，随机采样矩阵Ω′既包含对高频分量系数的随机放射线采样，又含有对大于阈值的低频系数的确定性采样。 (2) 将上述局部余弦变换系数作为观测采样值，通过求解式(4)在小波稀疏变换域中l₁范数最小化以及时域全变差正则化的问题重构原始地磁数据。

实现数据分频压缩与重构方法的部分伪码如下：

步骤1：加载地磁异常数据x，对数据进行归一化处理，设置二维网格数据的间距为Δx、Δy，网格数为m、n参数。

步骤2：进行余弦变换系数分频。对整体地磁数据进行离散余弦变换u=DCT(x)，设置分频系数阈值参数α，将变换域系数分解为高频和低频系数。

步骤3：局部余弦变换系数采样压缩。对低频系数进行无损压缩编码u′_L=u_L，利用放射线采样矩阵R对高频系数进行稀疏采样压缩编码u′H=R·u_H。

步骤4：获得重构的线性观测采样数据。对压缩编码系数解码，根据采样位置生成局部余弦变换测量观测值y=C_Ω·x。

步骤5：公式(4)可转换为拉格朗日无约束项形式的TVL1-L2压缩感知模型公式

点击浏览公式

步骤6：根据公式(6)，利用文献^[25]提出的RecPF算法求解TV正则化和l₁范数最小化的优化问题，最终获得压缩重构的地磁数据。该方法基于经典的带参数拉格朗日算法，利用方向选择算法实现优化问题的快速计算，在每次迭代过程中仅包含局部余弦变换系数阈值计算，迭代次数更少，抗差性更强。

通过步骤1～3实现了分频压缩，对大多数地磁数据而言，高频分量可近似认为是稀疏的，但低频分量不是稀疏的，它是数据主要能量在较低频率下的逼近信号。如果将高、低频系数一起与放射线采样矩阵R相乘，则在较低采样率的情况下会破坏低频分量系数之间的相关性，导致重构效果变差。另外由于低频系数所占的比重很小，约3%左右，因此完整保留这些少量的含有重要信息的系数不会影响到数据的压缩效率，反而由于精确地保留了这些系数，可以有效提高重构数据的质量。

而采用步骤4～6利用压缩感知方法精确重构地磁数据。该方法不仅在小波变换域中严格稀疏，而且在时域具有较小的全变差，有效地解决了数据平滑与轮廓维持之间的矛盾，在优化重构过程中抑制高频误差噪声的同时不对最优解加强平滑作用，较好地保留了地磁数据的轮廓和纹理特征，减少了重构中的高频误差。算法结构框图如图 1所示。

图 1 基于压缩感知的数据分频压缩及重构算法结构框图 Fig. 1 The block diagram of data frequency divisive compression and reconstruction algorithm based on compressed sensing

图选项

美国国家地球物理数据中心(NGDC)网站提供反映当地区域内磁场高频分量的地磁异常磁数据。以2002年发布的美国加州地区某处航磁测量数据作为本文仿真试验数据，机载磁力仪精度为1 nT，数据采样频率为2 Hz，地磁异常数据已消除了地磁场的日变影响，误差在2 nT以下或更低。对航磁测线数据经过Kriging插值、滤波处理形成了100像素×100像素网格数据，如图 2(a)所示，数据经度范围：-122.04°～-121.9°，纬度范围：36.20°～36.34°，分辨率为5″，地磁异常数据统计特征如表 1中所示。

图 2 原始地磁异常数据及余弦变换局部系数采样模式 Fig. 2 The original geomagnetic anomaly data and DCT domain sampling pattern

图选项

为分析重构数据精度，以及压缩感知和分频压缩在数据重构中独立发挥的作用，另外给出两个压缩重构方法进行对比：一个为滤波反投影方法(filtered back projection)；另一个为直接放射线采样方法，即对余弦变换系数不进行分频采样压缩，并同样利用压缩感知方法来重构地磁异常数据。

首先对地磁异常数据归一化处理并进行离散余弦变换，由于离散余弦变换具有频率分解和能量压缩的特性，使得数据在变换域中较大的系数集中在图像的左上角，本文方法通过设置离散余弦变换阈值参数α，将地磁数据在频域中分为低频和高频分量。本仿真试验中根据实际情况取α=0.5，则余弦变换低频分量系数采样率为2.45%，如图 2(b)所示。然后对高频分量系数进行随机采样，采样系数为28条放射线测量，如图 2(c)所示。将以上获得分频压缩总采样率为25%的采样系数作为本文方法及滤波反投影方法的重构观测数据。同时利用放射线直接采样获得同样总采样率的系数作为直接放射线采样方法重构观测数据。表 2给出不同采样方法高频和低频系数所占比例。

利用图像质量客观评价标准：峰值信噪比和相对误差对复原数据结果进行评价。峰值信噪比定义如下

点击浏览公式

点击浏览公式

图 3中(a)～(f)为总采样率均为25%时，3种压缩重构方法恢复地磁数据以及相对于原始数据误差。表 2给出3种重构方法恢复数据的峰值信噪比和相对误差统计结果。可以看出，虽然滤波反投影重构方法与分频压缩方法为相同的观测重构系数，但在重构算法中利用压缩感知能更精确地重构高频分量数据，峰值信噪比提高了4.81 dB。而与放射线采样方法相比，分频压缩方法通过阈值采样，在较低采样率下保留更多的低频分量系数，重构数据中含有更精确的低频分量数据信息，峰值信噪比相对放射线采样方法提高了2.97 dB。从重构算法的计算时间来看，本文方法需要求解l₁范数最小化和TV正则化等凸优化问题，虽然计算时间比滤波反投影重构方法长，但在可允许的计算时间内提高重构数据精度是值得的。

图 3 采样率为25%时3种重构方法恢复数据结果 Fig. 3 The retrieval results of three different algorithms with sampling ratio 25%

图选项

离散余弦变换阈值参数α的设置决定了在数据压缩中低频分量的采样率，而放射线采样矩阵R决定高频分量的采样率，为分析不同采样率对3种方法重构数据的峰值信噪比影响，图 4给出各方法在不同采样率下的峰值信噪比，以及相对于滤波反投影重构方法，其他两种方法增加的峰值信噪比。通过分析可以得出以下两点结论：

图 4 不同数据采样率下重构算法峰值信噪比分析 Fig. 4 Analysis of reconstruction algorithms’ PSNR via different sampling ratios

图选项

(1) 基于压缩感知的分频压缩方法在不同采样率下重构数据均具有更高的峰值信噪比，与传统滤波反投影方法相比可分别增加1～6 dB左右。表明将压缩感知方法应用于高压缩比的地磁数据压缩重构中具有明显的优势。

(2) 由于采用余弦变换阈值进行分频压缩，有效保留了低频分量数据，可以避免如放射线采样方法在采样率小于23%时，峰值信噪比低于滤波反投影重构方法的情况，表明本文方法在较低采样率条件下仍能保证较高的重构精度。随着采样率增加，放射线采样方法逐渐包含更多重要的低频分量系数，重构精度增加，在采样率大于40%以后精度与本文方法相同。

地磁辅助导航系统中需存储大量高精度地磁数据作为基准图数据库。针对传统数据压缩方法存在高频细节信息损失严重的问题，提出基于压缩感知理论的分频压缩联合重构方法，在余弦变换压缩编码过程中不对其低频变换系数压缩使之处于特殊的地位，而对高频变换系数进行稀疏采样压缩编码，实现数据压缩。这样不仅能较完整地保存低频分量数据，同时利用压缩感知重构算法，优化了重构方法中的随机观测矩阵，可以精确重构原始数据中高频分量特征，因此重构数据不会像滤波方法使整个数据图像平滑，而是保留清晰的轮廓，从而在较高数据压缩比条件下获得更高质量的重构数据。