基于Curvelet變換冗余字典的重力數據稀疏表示與重建

牛麗琨， 吳美平

(國防科學技術大學 機電工程與自動化學院，長沙 410073)

0 引言

全球重力數據作為各領域學科的基礎性資料存在很大空白，壓縮感知理論為實質上是欠奈奎斯特采樣的重力數據測量重構提供了理論指導[2]。由觀測量通過求解最優化問題得到重力數據的稀疏表示，實現由于客觀條件限制而極不完整、規則的重力測量數據重構[6]。該理論以信號的稀疏表示為前提，變換越稀疏，觀測矩陣與稀疏字典的相關性越小，重構的數據越精確[5]。

楊亞鵬[1]首次將壓縮感知理論應用于重力測量數據重構，選用了快速傅里葉變換基、離散小波基研究了重力數據的稀疏性。這些都屬于完備正交基展開法。傅里葉變換適用于一維重力數據，小波變換適用于二維重力數據，但只具備水平、垂直、對角線三個方向的選擇性，在處理復雜的重力場時有很大的局限性，不利于邊緣和紋理的捕捉處理。

構造合適的變換基，實現二維重力數據更為稀疏的表示對于測量矩陣的設計和重構精度的提高有很大意義[7-8]。筆者采用冗余字典作為變換基[3-4]，鑒于Curvelet變換具有非常好的各向異性特點，運用Curvelet變換作為核函數構造過完備冗余字典[9]。基于EGM2008全球重力場模型，對某海域的重力數據通過正交匹配算法進行稀疏分解，最后將得到的稀疏系數重建重力數據。

1 基于Curvelet變換得到冗余字典

運用Curvelet變換作為核函數構造過完備冗余字典。第一代Cuevelet變換實際上是子帶分解和脊波變換，考慮到計算機實現的快速和便利，筆者采用Candes、Donoho[2-4]提出的第二代Cuevelet變換中基于Wrap(Wrapping-based Transform)的快速離散曲波變換實現方法。

1.1 Wrap算法

Wrap算法的核心思想是圍繞原點Wrap，是指在具體實現時對任意區域，通過周期化技術一一映射到原點的仿射區域。運算速度快，算法效率高，具體過程如下：

1)對于給定的一個笛卡爾坐標系下的二維函數f進行二維快速傅里葉變換，得到二維頻域表示：

(1)

(2)

(3)

4)圍繞原點Wrapping局部化式(3)可得：

(4)

Curvelet變換需要對尺度和方向參數離散化，尺度數的劃分依據數據規模有所不同，將二維離散傅里葉平面劃分為一系列的同心矩形，每一個矩形環代表一個尺度。然后進行角度劃分，代表不同方向的曲波，如圖1所示，陰影部分即第四尺度層上第九個角度上的曲波。

圖1 尺度角度劃分Fig.1 Division of scale and angle

1.2 冗余字典

采用EGM2008全球重力場模型，計算了太平洋海域上N11°～N13°,E°133°～E°135°，分辨率為1′×1′的重力異常數據。如圖2所示，針對此121×121的重力圖在Curvelet冗余字典上進行稀疏表示。

圖2 某海域重力等值線圖Fig.2 Gravity contour map of a certain sea

此海域重力數據分布經過Curvelet變換得到C{j}{l}(k1,k2)結構的系數，j表示尺度，l表示方向，(k1,k2)表示尺度層上第l個方向的矩陣坐標，反映了空間位置信息。該重力圖經過Wrap算法后的曲波系數格式如表1所示。

121×121的重力數據將尺度劃分為4層，按照頻率由低到高為最內層的Coarse層(粗尺度層)，是低頻系數，包含了圖像概貌。中間Detail層(細尺度層)，包含邊緣特征。最外層的Fine層(精細尺度層)，是高頻系數，包含了細節、邊緣特征。

表1 Curvelet系數結構

不同尺度層不同的方向參數構成了冗余字典的所有原子，構造冗余字典時，選擇Coarse層和Fine層系數作為每部字典都有的原子，Detail層兩種尺度不同角度的曲波系數作為構成不同字典的不同原子。根據Detail層尺度和方向的不同，該重力圖共可以得到16+32=48種字典。

不同尺度下，不同方向的曲波系數不同，有的幅值較大，有的接近于“0”。分析后可以得到如下結論：①同一尺度下，方向與圖像邊緣相同的曲波系數較大；②同一方向上，尺度與捕捉對象相當的曲波在逼近中有更大的貢獻；③不同尺度、不同角度的曲波疊加起來構成完整的圖像。

圖3 Curvelet系數顯示(j=3,l=5 )Fig.3 Display of Curvelet coefficient when j=3,l=5

對頻率域的曲波作逆變換得到空間域的波形，圖3給出第3尺度第5個角度的笛卡爾坐標下Curvelet系數三維顯示。

隨著尺度和角度的不同，在這48個原子中，字典重建的效果會有差異。比較這些原子，可以得到更好的重建效果。用n表示選取的最優原子個數。

1.3 稀疏表示

使用正交匹配追蹤(OMP)算法，對重力數據與冗余字典進行OMP計算，得到重力數據在Curvelet冗余字典下最佳的稀疏表示。步驟為：

1)令初始殘差信號為元信號，初始迭代次數t=0，根據計算得到的最優原子個數n和所需稀疏分解的精確度確定最大原子個數m。

2)計算殘差信號與冗余字典中所有原子的內積，并從冗余字典中，選擇使得內積最大的原子。

3)利用施密特正交化方法對內積最大原子的乘積來更新殘差信號。

4)判斷是否迭代，如果次數t>m，則停止，否則繼續迭代。

2 圖像的重建和效果評估

2.1 圖像稀疏重建

用不同原子構成的冗余字典去重建原始重力數據，與原重力數據作對比。這里選擇第2尺度下第2個角度的原子構成的字典重建重力數據，結果如圖4所示。

圖4 重建重力圖(j=2,l=2)Fig.4 Gravity map of reconstruction when j=2,l=2

2.2 效果評估

為了比較這些原子，對重建的數據設立指標進行評估。

2.2.1 PSNR值

峰值信噪比(PSNR)，一種評價圖像的客觀標準，單位是dB，計算公式為式(5)。

(5)

其中：MSE為原重力數據與重建重力數據之間的均方誤差；n是每個采樣值的比特數。PSNR值越大，就代表重建效果越好。所以計算重建圖像的PSNR，作為衡量重建圖像的標準之一。

2.2.2 稀疏度φ值

稀疏度φ就是稀疏系數矩陣非零項個數占原始大小的比例，稀疏度越低，圖像表示越稀疏。

2.2.3η值

計算第2尺度下第2個角度的原子構成的字典重建效果各評價標準參數。PSNR=43.821 9 dB，稀疏度φ=59.51%。

3 基于Curvelet變換的字典優化

對原始重力數據，用不同原子組成的字典進行重建，根據不同尺度不同方向的原子測試重建效果。

如圖5所示，為48種不同原子重建重力數據的PSNR值和稀疏度φ值。

根據PSNR值的大小作為標準選擇最優原子時，最優原子是來自第6、7、14、15、22、23、30、31個原子，得到表2，重力數據重建后的PSNR值和稀疏度的值。

根據表2，對最優原子進行組合，得到新的字典，對原始重力數據進行重建，得到PSNR值和稀疏度,如表3所示。

圖5 48種原子組成Curvelet字典重建結果Fig.5 Reconstruction of 48 kinds Curvelet dictionary(a) PSNR值；(b)φ值

第n項67141522233031PSNR/dB44.011144.010744.011144.010744.011144.010738.3638.36稀疏度/%37.6037.6037.6037.6038.3638.3660.7360.73

表3 最優原子組合重建效果

從表3可以看出，隨著最優原子的增加，PSNR逐步增大，稀疏度也逐漸增大。

采用另一個標準，即η=PSNR/φ來判斷，PSNR越大，稀疏度越小，η就越大。計算48種不同原子重建重力數據的η值(圖6)。根據圖6可以看出，η最大的是第2、3、6、7、10、11、14、15個原子，得到表4，重力數據重建后的η值。

根據表4，對最優原子進行組合，得到新的字典，對原始重力數據進行重建，得到PSNR值、稀疏度和η值,如表5所示。

圖6 48種原子組成Curvelet字典重建η變化Fig.6 change of η by 48 kinds Curvelet dictionary

第n項236710111415η值1.1651.1641.1711.171.16561.16391.17061.1706

表5 最優原子組合重建效果

將表5的結果反映在圖像上，如圖7所示。

由圖7可以看出，當原子組合數是3的時候，獲得的η值最高。因此原子組合在n=3的時候最佳。

圖7 最優原子組合重建效果曲線圖Fig.7 The graph of optimal atomic combination

圖8 重建重力數據圖Fig.8 The graph of gravity reconstruction

參考最優原子個數，選擇5作為最大原子個數，按照OMP算法對重力數據進行稀疏表示。最終稀疏表示所選原子為：2、6、8、14、15。稀疏度為42.98%，重建重力數據如圖8所示。

圖9 重力數據對比(E=133.5°)Fig.9 The contrast of gravity data when E=133.5°

圖10 重力數據對比(E=134°)Fig.10 The contrast of gravity data when E=134°

為了更加直觀地反映重建的精度效果，將圖2和圖8沿經線做剖面對比。文獻[1]中采用小波基稀作為稀疏變換基，同樣使用OMP算法進行重建，針對不同的重力區域，稀疏度均在50%以上，重建重力圖像峰值信噪比(PSNR)值在30左右，其η值在0.6左右。

由圖9、圖10可以看出對于二維重力數據，相比于文獻[1]，我們所采用的壓縮感知稀疏表示與重建算法很大程度上降低了數據的稀疏度和重建的精度，對于使用壓縮感知算法進行重力數據測量具有重大意義。

4 結論

筆者利用過完備冗余字典作為變換基，實現了二維重力數據相對小波基更為稀疏的信號表示，并用得到的稀疏系數重建了重力數據。

鑒于Curvelet變換具有非常好的各向異性特點，采用Curvelet變換作為核函數構造冗余字典。選取不同尺度不同方向的曲波系數作為原子，構造不同的字典對重力數據進行稀疏表示和重建，并對這些原子進行比較，用PSNR值、稀疏度和η值作為評估標準確定最優原子個數。最后采用正交匹配算法對重力數據進行稀疏表示，并用得到的稀疏系數重建了重力數據。