劉繁明，張迎發(fā)，李 艷，姚建奇

哈爾濱工程大學 自動化學院，黑龍江 哈爾濱 150001

1 引 言

隨著航空地磁測量精度的提高以及高精度插值算法的應用，地磁數(shù)據(jù)的網(wǎng)格分辨率達到百米量級左右，大大提高了數(shù)據(jù)的空間分辨率，其豐富的紋理信息，更利于應用地磁及其梯度數(shù)據(jù)進行輔助導航［1－5］。但由于數(shù)據(jù)量的劇增，對導航數(shù)據(jù)庫構(gòu)建提出了較高的要求，研究大容量地磁數(shù)據(jù)的壓縮方法，降低高分辨率數(shù)據(jù)庫的存儲容量具有重要意義。

傳統(tǒng)數(shù)據(jù)壓縮方法是以犧牲數(shù)據(jù)的高頻細節(jié)信息為代價，通常運用BDCT（blocking discrete cosine transform）編解碼技術(shù)對數(shù)據(jù)進行分塊變換處理［6］。較大分塊數(shù)據(jù)相關(guān)性較高，可以實現(xiàn)較大的壓縮比，但大尺寸分塊數(shù)據(jù)壓縮時又會丟失數(shù)據(jù)的平穩(wěn)性，引入振鈴效應（ringing effect）噪聲和塊效應（block effect）邊界誤差［7－9］。由此重構(gòu)的數(shù)據(jù)高頻輪廓特征信息損失嚴重或引入重構(gòu)噪聲，增大了地磁輔助導航系統(tǒng)匹配失敗或是誤匹配的概率［10－11］，因此不適用于大容量的導航數(shù)據(jù)壓縮。

壓縮感知（compressed sensing，CS）作為一種新興的稀疏信號精確重構(gòu)理論已成為核磁共振和合成孔徑雷達成像等圖像處理領(lǐng)域的重要方法［12－13］。由于地磁數(shù)據(jù)能量多集中在低頻分量上，而壓縮感知對稀疏信號的隨機采樣矩陣通常為局部傅里葉變換系數(shù)星形采樣，這種采樣矩陣破壞了低頻分量系數(shù)之間的相關(guān)性，影響了數(shù)據(jù)重構(gòu)精度［14－15］。

為解決大容量地磁數(shù)據(jù)壓縮的同時，又能較好地保存細節(jié)輪廓信息的問題，本文提出一種基于壓縮感知理論的數(shù)據(jù)分頻壓縮及聯(lián)合重構(gòu)方法。

2 壓縮感知理論

壓縮感知是一種新興理論，它根據(jù)信號本身或在變換域中具有的稀疏特性，利用遠少于Nyquist采樣定理所需的采樣數(shù)據(jù)就可以精確地重構(gòu)原始數(shù)據(jù)，因此，該理論在各個領(lǐng)域中展開廣泛研究及應用［16］。

根據(jù)壓縮感知理論，當隨機測量矩陣Φ與稀疏變換矩陣Ψ滿足約束等距性（restricted isometry property，RIP）條件，即可利用信號x線性測量獲得的觀測值y＝Φx＝Au，通過解決以下l1范數(shù)最小值問題來實現(xiàn)原始數(shù)據(jù)的精確重構(gòu)

式中，A＝ΦΨ為感知矩陣，當線性測量值含有高斯噪聲時，根據(jù)噪聲大小σ對上式中l(wèi)1范數(shù)約束條件適當松弛，尋找l1范數(shù)最小值問題可以表述為以下線性凸優(yōu)化問題

對于時域中絕大多數(shù)圖像數(shù)據(jù)，其離散梯度具有稀疏性，全變差（total variation，TV）正則化模型能夠比較合理地保持圖像的輪廓結(jié)構(gòu)和邊緣特征［17］。因此，全變差方法常作為稀疏變換，替代l1范數(shù)應用在壓縮感知方法中［12］。文獻［18—19］結(jié)合TV正則化和l1范數(shù)的優(yōu)點，給出同時在稀疏變換域以及稀疏梯度時域中尋找最優(yōu)解的壓縮感知模型

式中，規(guī)劃參數(shù)λ≥0，選擇適當?shù)摩丝梢哉{(diào)整稀疏變換域中l(wèi)1范數(shù)和時域全變差兩種方法在凸優(yōu)化尋求最優(yōu)解過程中所占的比重。對于圖像數(shù)據(jù)x，其全變差定義為

式中，D1x＝x（i＋1，j）－x（i，j），D2x＝x（i，j＋1）－x（i，j）為圖像的像素點x（i，j）在水平和垂直方向上的有限差分。

由于TV正則化和l1范數(shù)都具有非平滑、不可微分的特性，傳統(tǒng)優(yōu)化算法如共軛梯度法和匹配追蹤算法在解決以上優(yōu)化問題時，需要耗費大量的計算 時 間［20－21］。 近 期 一 些 快 速 算 法 如 TwIST、FISTA、TVAL3和RecPF等的出現(xiàn)，使得在求解壓縮感知問題時可以節(jié)省大量的時間［22－25］。

3 基于壓縮感知的分頻壓縮及重構(gòu)方法

對于導航用地磁數(shù)據(jù)，壓縮時需要保留精度較高的邊緣、輪廓和紋理等局部高頻特征，基于壓縮感知理論的數(shù)據(jù)分頻壓縮及聯(lián)合重構(gòu)方法主要包含兩部分內(nèi)容：頻域余弦變換的分頻壓縮和基于壓縮感知的聯(lián)合重構(gòu)。

分頻壓縮方法的主要過程為：

（1）以較大范圍的高精度地磁數(shù)據(jù)作為整體進行二維離散余弦變換，以將數(shù)據(jù)進行多尺度頻率分解。

（2）設置變換域系數(shù)閾值α，將幅值大于閾值的低頻分量系數(shù)進行無損編碼。

（3）高頻分量系數(shù)包含邊緣和紋理成分，利用放射線采樣矩陣R對高頻分量系數(shù)進行稀疏采樣編碼。

數(shù)據(jù)重構(gòu)方法的主要過程為：

（1）對分頻壓縮存儲的低頻系數(shù)和高頻系數(shù)按所處位置進行重組，作為壓縮感知重構(gòu)中的線性觀測值y＝Φx。其中，測量矩陣Φ為局部余弦變換系數(shù)矩陣Φ＝CΩ＝Ω′·C［·］，C［·］代表二維離散余弦變換，隨機采樣矩陣Ω′既包含對高頻分量系數(shù)的隨機放射線采樣，又含有對大于閾值的低頻系數(shù)的確定性采樣。

（2）將上述局部余弦變換系數(shù)作為觀測采樣值，通過求解式（4）在小波稀疏變換域中l(wèi)1范數(shù)最小化以及時域全變差正則化的問題重構(gòu)原始地磁數(shù)據(jù)。

實現(xiàn)數(shù)據(jù)分頻壓縮與重構(gòu)方法的部分偽碼如下：

步驟1：加載地磁異常數(shù)據(jù)x，對數(shù)據(jù)進行歸一化處理，設置二維網(wǎng)格數(shù)據(jù)的間距為Δx、Δy，網(wǎng)格數(shù)為m、n參數(shù)。

步驟2：進行余弦變換系數(shù)分頻。對整體地磁數(shù)據(jù)進行離散余弦變換u＝DCT（x），設置分頻系數(shù)閾值參數(shù)α，將變換域系數(shù)分解為高頻和低頻系數(shù)。

步驟3：局部余弦變換系數(shù)采樣壓縮。對低頻系數(shù)進行無損壓縮編碼u′L＝uL，利用放射線采樣矩陣R對高頻系數(shù)進行稀疏采樣壓縮編碼u′H＝R·uH。

步驟4：獲得重構(gòu)的線性觀測采樣數(shù)據(jù)。對壓縮編碼系數(shù)解碼，根據(jù)采樣位置生成局部余弦變換測量觀測值y＝CΩ·x。

步驟5：公式（4）可轉(zhuǎn)換為拉格朗日無約束項形式的TVL1－L2壓縮感知模型公式

步驟6：根據(jù)公式（6），利用文獻［25］提出的RecPF算法求解TV正則化和l1范數(shù)最小化的優(yōu)化問題，最終獲得壓縮重構(gòu)的地磁數(shù)據(jù)。該方法基于經(jīng)典的帶參數(shù)拉格朗日算法，利用方向選擇算法實現(xiàn)優(yōu)化問題的快速計算，在每次迭代過程中僅包含局部余弦變換系數(shù)閾值計算，迭代次數(shù)更少，抗差性更強。

通過步驟1～3實現(xiàn)了分頻壓縮，對大多數(shù)地磁數(shù)據(jù)而言，高頻分量可近似認為是稀疏的，但低頻分量不是稀疏的，它是數(shù)據(jù)主要能量在較低頻率下的逼近信號。如果將高、低頻系數(shù)一起與放射線采樣矩陣R相乘，則在較低采樣率的情況下會破壞低頻分量系數(shù)之間的相關(guān)性，導致重構(gòu)效果變差。另外由于低頻系數(shù)所占的比重很小，約3%左右，因此完整保留這些少量的含有重要信息的系數(shù)不會影響到數(shù)據(jù)的壓縮效率，反而由于精確地保留了這些系數(shù)，可以有效提高重構(gòu)數(shù)據(jù)的質(zhì)量。

而采用步驟4～6利用壓縮感知方法精確重構(gòu)地磁數(shù)據(jù)。該方法不僅在小波變換域中嚴格稀疏，而且在時域具有較小的全變差，有效地解決了數(shù)據(jù)平滑與輪廓維持之間的矛盾，在優(yōu)化重構(gòu)過程中抑制高頻誤差噪聲的同時不對最優(yōu)解加強平滑作用，較好地保留了地磁數(shù)據(jù)的輪廓和紋理特征，減少了重構(gòu)中的高頻誤差。算法結(jié)構(gòu)框圖如圖1所示。

圖1 基于壓縮感知的數(shù)據(jù)分頻壓縮及重構(gòu)算法結(jié)構(gòu)框圖Fig.1 The block diagram of data frequency divisive compression and reconstruction algorithm based on compressed sensing

4 實測數(shù)據(jù)仿真試驗分析

美國國家地球物理數(shù)據(jù)中心（NGDC）網(wǎng)站提供反映當?shù)貐^(qū)域內(nèi)磁場高頻分量的地磁異常磁數(shù)據(jù)。以2002年發(fā)布的美國加州地區(qū)某處航磁測量數(shù)據(jù)作為本文仿真試驗數(shù)據(jù)，機載磁力儀精度為1nT，數(shù)據(jù)采樣頻率為2Hz，地磁異常數(shù)據(jù)已消除了地磁場的日變影響，誤差在2nT以下或更低。對航磁測線數(shù)據(jù)經(jīng)過Kriging插值、濾波處理形成了100像素×100像素網(wǎng)格數(shù)據(jù)，如圖2（a）所示，數(shù)據(jù)經(jīng)度范圍：－122.04°～－121.9°，緯度范圍：36.20°～36.34°，分辨率為5″，地磁異常數(shù)據(jù)統(tǒng)計特征如表1中所示。

表1 地磁異常數(shù)據(jù)統(tǒng)計特征Tab.1 Statistical character of the geomagnetic anomaly data nT

為分析重構(gòu)數(shù)據(jù)精度，以及壓縮感知和分頻壓縮在數(shù)據(jù)重構(gòu)中獨立發(fā)揮的作用，另外給出兩個壓縮重構(gòu)方法進行對比：一個為濾波反投影方法（filtered back projection）；另一個為直接放射線采樣方法，即對余弦變換系數(shù)不進行分頻采樣壓縮，并同樣利用壓縮感知方法來重構(gòu)地磁異常數(shù)據(jù)。

圖2 原始地磁異常數(shù)據(jù)及余弦變換局部系數(shù)采樣模式Fig.2 The original geomagnetic anomaly data and DCT domain sampling pattern

首先對地磁異常數(shù)據(jù)歸一化處理并進行離散余弦變換，由于離散余弦變換具有頻率分解和能量壓縮的特性，使得數(shù)據(jù)在變換域中較大的系數(shù)集中在圖像的左上角，本文方法通過設置離散余弦變換閾值參數(shù)α，將地磁數(shù)據(jù)在頻域中分為低頻和高頻分量。本仿真試驗中根據(jù)實際情況取α＝0.5，則余弦變換低頻分量系數(shù)采樣率為2.45%，如圖2（b）所示。然后對高頻分量系數(shù)進行隨機采樣，采樣系數(shù)為28條放射線測量，如圖2（c）所示。將以上獲得分頻壓縮總采樣率為25%的采樣系數(shù)作為本文方法及濾波反投影方法的重構(gòu)觀測數(shù)據(jù)。同時利用放射線直接采樣獲得同樣總采樣率的系數(shù)作為直接放射線采樣方法重構(gòu)觀測數(shù)據(jù)。表2給出不同采樣方法高頻和低頻系數(shù)所占比例。

利用圖像質(zhì)量客觀評價標準：峰值信噪比和相對誤差對復原數(shù)據(jù)結(jié)果進行評價。峰值信噪比定義如下

相對誤差定義為

式中，fmax（x，y）為地磁異常數(shù)據(jù)最大值；f（x，y）為原始數(shù)據(jù)；f′（x，y）為待評價重構(gòu)數(shù)據(jù)；M、N為數(shù)據(jù)網(wǎng)格數(shù)；‖·‖2為l2范數(shù)。

圖3 采樣率為25%時3種重構(gòu)方法恢復數(shù)據(jù)結(jié)果Fig.3 The retrieval results of three different algorithms with sampling ratio 25%

圖3中（a）～（f）為總采樣率均為25%時，3種壓縮重構(gòu)方法恢復地磁數(shù)據(jù)以及相對于原始數(shù)據(jù)誤差。表2給出3種重構(gòu)方法恢復數(shù)據(jù)的峰值信噪比和相對誤差統(tǒng)計結(jié)果。可以看出，雖然濾波反投影重構(gòu)方法與分頻壓縮方法為相同的觀測重構(gòu)系數(shù)，但在重構(gòu)算法中利用壓縮感知能更精確地重構(gòu)高頻分量數(shù)據(jù)，峰值信噪比提高了4.81dB。而與放射線采樣方法相比，分頻壓縮方法通過閾值采樣，在較低采樣率下保留更多的低頻分量系數(shù)，重構(gòu)數(shù)據(jù)中含有更精確的低頻分量數(shù)據(jù)信息，峰值信噪比相對放射線采樣方法提高了2.97dB。從重構(gòu)算法的計算時間來看，本文方法需要求解l1范數(shù)最小化和TV正則化等凸優(yōu)化問題，雖然計算時間比濾波反投影重構(gòu)方法長，但在可允許的計算時間內(nèi)提高重構(gòu)數(shù)據(jù)精度是值得的。

表2 重構(gòu)算法恢復效果及其誤差統(tǒng)計結(jié)果對比Tab.2 The comparison of three reconstruction algorithms and their error statistics results

離散余弦變換閾值參數(shù)α的設置決定了在數(shù)據(jù)壓縮中低頻分量的采樣率，而放射線采樣矩陣R決定高頻分量的采樣率，為分析不同采樣率對3種方法重構(gòu)數(shù)據(jù)的峰值信噪比影響，圖4給出各方法在不同采樣率下的峰值信噪比，以及相對于濾波反投影重構(gòu)方法，其他兩種方法增加的峰值信噪比。通過分析可以得出以下兩點結(jié)論：

（1）基于壓縮感知的分頻壓縮方法在不同采樣率下重構(gòu)數(shù)據(jù)均具有更高的峰值信噪比，與傳統(tǒng)濾波反投影方法相比可分別增加1～6dB左右。表明將壓縮感知方法應用于高壓縮比的地磁數(shù)據(jù)壓縮重構(gòu)中具有明顯的優(yōu)勢。

（2）由于采用余弦變換閾值進行分頻壓縮，有效保留了低頻分量數(shù)據(jù)，可以避免如放射線采樣方法在采樣率小于23%時，峰值信噪比低于濾波反投影重構(gòu)方法的情況，表明本文方法在較低采樣率條件下仍能保證較高的重構(gòu)精度。隨著采樣率增加，放射線采樣方法逐漸包含更多重要的低頻分量系數(shù)，重構(gòu)精度增加，在采樣率大于40%以后精度與本文方法相同。

圖4 不同數(shù)據(jù)采樣率下重構(gòu)算法峰值信噪比分析Fig.4 Analysis of reconstruction algorithms’PSNR via different sampling ratios

5 結(jié) 論

地磁輔助導航系統(tǒng)中需存儲大量高精度地磁數(shù)據(jù)作為基準圖數(shù)據(jù)庫。針對傳統(tǒng)數(shù)據(jù)壓縮方法存在高頻細節(jié)信息損失嚴重的問題，提出基于壓縮感知理論的分頻壓縮聯(lián)合重構(gòu)方法，在余弦變換壓縮編碼過程中不對其低頻變換系數(shù)壓縮使之處于特殊的地位，而對高頻變換系數(shù)進行稀疏采樣壓縮編碼，實現(xiàn)數(shù)據(jù)壓縮。這樣不僅能較完整地保存低頻分量數(shù)據(jù)，同時利用壓縮感知重構(gòu)算法，優(yōu)化了重構(gòu)方法中的隨機觀測矩陣，可以精確重構(gòu)原始數(shù)據(jù)中高頻分量特征，因此重構(gòu)數(shù)據(jù)不會像濾波方法使整個數(shù)據(jù)圖像平滑，而是保留清晰的輪廓，從而在較高數(shù)據(jù)壓縮比條件下獲得更高質(zhì)量的重構(gòu)數(shù)據(jù)。

［1］ HEMANT K，THéBAULT E，MANDEA M，et al.Magnetic Anomaly Map of the World：Merging Satellite，Airborne，Marine and Ground－based Magnetic Data Sets［J］.Earth and Planetary Science Letters，2007，260（1－2）：56－71.

［2］ NABIGHIAN M N，GRAUCH V J S，HANSEN R O，et al.The Historical Development of the Magnetic Method in Exploration［J］.Geophysics，2005，70（6）：33N－61N.

［3］ ZHOU Jun，GE Zhilei，SHI Guiguo，et al.Key Technique and Development for Geomagnetic Navigation［J］.Journal of Astronautics，2008，29（5）：1467－1472.（周軍，葛致磊，施桂國，等.地磁導航發(fā)展與關(guān)鍵技術(shù)［J］.宇航學報，2008，29（5）：1467－1472.）

［4］ QIAN Dong，LIU Fanming，LI Yan，et al.Comparison of Gravity Gradient Reference Map Composition for Navigation［J］.Acta Geodaetica et Cartographica Sinica，2011，40（6）：736－744.（錢東，劉繁明，李艷，等.導航用重力梯度基準圖構(gòu)建方法的比較研究［J］.測繪學報，2011，40（6）：736－744.）

［5］ WANG Shicheng，WANG Zhe，ZHANG Jinsheng，et al.Technology of Preparation of Reference Map Using Total Geomagnetic Intensity Gradient Module［J］.Systems Engineering and Electronics，2009，31（4）：881－885.（王仕成，王哲，張金生，等.總磁場強度梯度模作為匹配特征量的基準圖制備技術(shù)［J］.系統(tǒng)工程與電子技術(shù)，2009，31（4）：881－885.）

［6］ RAO K R，YIP P.Discrete Cosine Transform：Algorithms，Advantages， Applications ［M］.Waltham：Academic Press，1990.

［7］ KASEZAWA T.Blocking Artifacts Reduction Using Discrete Cosine Transform［J］.IEEE Transactions on Consumer Electronics，1997，43（1）：48－55.

［8］ XU Jinbo.An Improved Algorithm to Remove the Block Effect in Transform Encoder Images［J］.Computer Engineering and Science，2006，28（2）：51－53.（徐金波.變換編碼中消除圖像“塊效應”的優(yōu)化算法［J］.計算機工程與科學，2006，28（2）：51－53.）

［9］ ZOU J J，YAN H.A Deblocking Method for BDCT Compressed Images Based on Adaptive Projections［J］.IEEE Transactions on Circuits and Systems for Video Technology，2005，15（3）：430－435.

［10］ LUO S，WANG Y，LIU Y，et al.Geomagnetic－Matching Technology Based on Improved ICP Algorithm［J］.Journal of Computer Applications，2008（s1）：351－354.

［11］ WANG Xianglei.The Study of Some Key Technology in Geomagnetic Matching Navigation［J］.Engineering of Surveying and Mapping，2011，20（1）：1－5.（王向磊.地磁匹配導航中幾項關(guān)鍵技術(shù)研究［J］.測繪工程，2011，20（1）：1－5.）

［12］ CANDéS 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.

［13］ PATEL V M，EASLEY G R，HEALY J D M，et al.Compressed Synthetic Aperture Radar［J］.IEEE Journal on Selected Topics in Signal Processing，2010，4（2）：244－254.

［14］ CANDES E J，TAO T.Near－optimal Signal Recovery from Random Projections：Universal Encoding Strategies？［J］.IEEE Transactions on Information Theory，2006，52（12）：5406－5425.

［15］ PATEL V M，MALEH R，GILBERT A C，et al.Gradientbased Image Recovery Methods from Incomplete Fourier Measurements［J］.IEEE Transactions on Image Processing，2012，21（1）：94－105.

［16］ TROPP J A，LASKA J N，DUARTE M F，et al.Beyond Nyquist：Efficient Sampling of Sparse Bandlimited Signals［J］.IEEE Transactions on Information Theory，2010，56（1）：520－544.

［17］ RUDIN L I，OSHER S，F(xiàn)ATEMI E.Nonlinear Total Variation Based Noise Removal Algorithms［J］.Physica D：Nonlinear Phenome－na，1992，60（1－4）：259－268.

［18］ CANDéS E J，ROMBERG J K.Signal Recovery from Random Projections［C］∥Proceedings of SPIE.San Jose：［s.n.］，2005.

［19］ CANDéS E，BRAUN N，WAKIN M.Sparse Signal and Image Recovery from Compressive Samples［C］∥Proceedings of International Symposium on Biomedical Imaging： From Nano to Macro.Arlington： ［s.n.］，2007.

［20］ LUSTIG M，DONOHO D，PAULY J M.Sparse MRI：The Application of Compressed Sensing for Rapid MR Imaging［J］.Magnetic Resonance in Medicine，2007，58（6）：1182－1195.

［21］ TROPP J A.Greed is Good：Algorithmic Results for Sparse Approximation［J］.IEEE Transactions on Information Theory，2004，50（10）：2231－2242.

［22］ BIOUCAS－DIAS J M，F(xiàn)IGUEIREDO M A T.A New TwIST： Two－step Iterative Shrinkage／Thresholding Algorithms for Image Restoration［J］.IEEE Transactions on Image Processing，2007，16（12）：2992－3004.

［23］ BECK A，TEBOULLE M.Fast Gradient－based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems［J］.IEEE Transactions on Image Processing，2009，18（11）：2419－2434.

［24］ ZUO W，LIN Z.A Generalized Accelerated Proximal Gradient Approach for Total－Variation－Based Image Restoration［J］.IEEE Transactions on Image Processing，2011，20（10）：2748－2759.

［25］ YANG J，ZHANG Y，YIN W.A Fast Alternating Direction Method for TVL1－L2Signal Reconstruction from Partial Fourier Data［J］.IEEE Journal on Selected Topics in Signal Processing，2010，4（2）：288－297.