基于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值、稀疏度和η值作为评估标准确定最优原子个数。最后采用正交匹配算法对重力数据进行稀疏表示，并用得到的稀疏系数重建了重力数据。