牛麗琨， 吳美平

(國(guó)防科學(xué)技術(shù)大學(xué) 機(jī)電工程與自動(dòng)化學(xué)院，長(zhǎng)沙 410073)

0 引言

全球重力數(shù)據(jù)作為各領(lǐng)域?qū)W科的基礎(chǔ)性資料存在很大空白，壓縮感知理論為實(shí)質(zhì)上是欠奈奎斯特采樣的重力數(shù)據(jù)測(cè)量重構(gòu)提供了理論指導(dǎo)[2]。由觀測(cè)量通過求解最優(yōu)化問題得到重力數(shù)據(jù)的稀疏表示，實(shí)現(xiàn)由于客觀條件限制而極不完整、規(guī)則的重力測(cè)量數(shù)據(jù)重構(gòu)[6]。該理論以信號(hào)的稀疏表示為前提，變換越稀疏，觀測(cè)矩陣與稀疏字典的相關(guān)性越小，重構(gòu)的數(shù)據(jù)越精確[5]。

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

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

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

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

1.1 Wrap算法

Wrap算法的核心思想是圍繞原點(diǎn)Wrap，是指在具體實(shí)現(xiàn)時(shí)對(duì)任意區(qū)域，通過周期化技術(shù)一一映射到原點(diǎn)的仿射區(qū)域。運(yùn)算速度快，算法效率高，具體過程如下：

1)對(duì)于給定的一個(gè)笛卡爾坐標(biāo)系下的二維函數(shù)f進(jìn)行二維快速傅里葉變換，得到二維頻域表示：

(1)

(2)

(3)

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

(4)

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

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

1.2 冗余字典

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

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

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

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

表1 Curvelet系數(shù)結(jié)構(gòu)

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

不同尺度下，不同方向的曲波系數(shù)不同，有的幅值較大，有的接近于“0”。分析后可以得到如下結(jié)論：①同一尺度下，方向與圖像邊緣相同的曲波系數(shù)較大；②同一方向上，尺度與捕捉對(duì)象相當(dāng)?shù)那ㄔ诒平杏懈蟮呢暙I(xiàn)；③不同尺度、不同角度的曲波疊加起來構(gòu)成完整的圖像。

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

對(duì)頻率域的曲波作逆變換得到空間域的波形，圖3給出第3尺度第5個(gè)角度的笛卡爾坐標(biāo)下Curvelet系數(shù)三維顯示。

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

1.3 稀疏表示

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

1)令初始?xì)埐钚盘?hào)為元信號(hào)，初始迭代次數(shù)t=0，根據(jù)計(jì)算得到的最優(yōu)原子個(gè)數(shù)n和所需稀疏分解的精確度確定最大原子個(gè)數(shù)m。

2)計(jì)算殘差信號(hào)與冗余字典中所有原子的內(nèi)積，并從冗余字典中，選擇使得內(nèi)積最大的原子。

3)利用施密特正交化方法對(duì)內(nèi)積最大原子的乘積來更新殘差信號(hào)。

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

2 圖像的重建和效果評(píng)估

2.1 圖像稀疏重建

用不同原子構(gòu)成的冗余字典去重建原始重力數(shù)據(jù)，與原重力數(shù)據(jù)作對(duì)比。這里選擇第2尺度下第2個(gè)角度的原子構(gòu)成的字典重建重力數(shù)據(jù)，結(jié)果如圖4所示。

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

2.2 效果評(píng)估

為了比較這些原子，對(duì)重建的數(shù)據(jù)設(shè)立指標(biāo)進(jìn)行評(píng)估。

2.2.1 PSNR值

峰值信噪比(PSNR)，一種評(píng)價(jià)圖像的客觀標(biāo)準(zhǔn)，單位是dB，計(jì)算公式為式(5)。

(5)

其中：MSE為原重力數(shù)據(jù)與重建重力數(shù)據(jù)之間的均方誤差；n是每個(gè)采樣值的比特?cái)?shù)。PSNR值越大，就代表重建效果越好。所以計(jì)算重建圖像的PSNR，作為衡量重建圖像的標(biāo)準(zhǔn)之一。

2.2.2 稀疏度φ值

稀疏度φ就是稀疏系數(shù)矩陣非零項(xiàng)個(gè)數(shù)占原始大小的比例，稀疏度越低，圖像表示越稀疏。

2.2.3η值

計(jì)算第2尺度下第2個(gè)角度的原子構(gòu)成的字典重建效果各評(píng)價(jià)標(biāo)準(zhǔn)參數(shù)。PSNR=43.821 9 dB，稀疏度φ=59.51%。

3 基于Curvelet變換的字典優(yōu)化

對(duì)原始重力數(shù)據(jù)，用不同原子組成的字典進(jìn)行重建，根據(jù)不同尺度不同方向的原子測(cè)試重建效果。

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

根據(jù)PSNR值的大小作為標(biāo)準(zhǔn)選擇最優(yōu)原子時(shí)，最優(yōu)原子是來自第6、7、14、15、22、23、30、31個(gè)原子，得到表2，重力數(shù)據(jù)重建后的PSNR值和稀疏度的值。

根據(jù)表2，對(duì)最優(yōu)原子進(jìn)行組合，得到新的字典，對(duì)原始重力數(shù)據(jù)進(jìn)行重建，得到PSNR值和稀疏度,如表3所示。

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

第n項(xiàng)67141522233031PSNR/dB44.011144.010744.011144.010744.011144.010738.3638.36稀疏度/%37.6037.6037.6037.6038.3638.3660.7360.73

表3 最優(yōu)原子組合重建效果

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

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

根據(jù)表4，對(duì)最優(yōu)原子進(jìn)行組合，得到新的字典，對(duì)原始重力數(shù)據(jù)進(jìn)行重建，得到PSNR值、稀疏度和η值,如表5所示。

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

第n項(xiàng)236710111415η值1.1651.1641.1711.171.16561.16391.17061.1706

表5 最優(yōu)原子組合重建效果

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

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

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

圖8 重建重力數(shù)據(jù)圖Fig.8 The graph of gravity reconstruction

參考最優(yōu)原子個(gè)數(shù)，選擇5作為最大原子個(gè)數(shù)，按照OMP算法對(duì)重力數(shù)據(jù)進(jìn)行稀疏表示。最終稀疏表示所選原子為：2、6、8、14、15。稀疏度為42.98%，重建重力數(shù)據(jù)如圖8所示。

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

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

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

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

4 結(jié)論

筆者利用過完備冗余字典作為變換基，實(shí)現(xiàn)了二維重力數(shù)據(jù)相對(duì)小波基更為稀疏的信號(hào)表示，并用得到的稀疏系數(shù)重建了重力數(shù)據(jù)。

鑒于Curvelet變換具有非常好的各向異性特點(diǎn)，采用Curvelet變換作為核函數(shù)構(gòu)造冗余字典。選取不同尺度不同方向的曲波系數(shù)作為原子，構(gòu)造不同的字典對(duì)重力數(shù)據(jù)進(jìn)行稀疏表示和重建，并對(duì)這些原子進(jìn)行比較，用PSNR值、稀疏度和η值作為評(píng)估標(biāo)準(zhǔn)確定最優(yōu)原子個(gè)數(shù)。最后采用正交匹配算法對(duì)重力數(shù)據(jù)進(jìn)行稀疏表示，并用得到的稀疏系數(shù)重建了重力數(shù)據(jù)。