CN102855644B - 一种基于小波自反馈的大规模地形数据压缩和解压方法 - Google Patents

Abstract

本发明针对大规模地形数据量大、应用背景复杂的问题，提出一种基于小波自反馈的大规模地形数据压缩和解压方法，该方法根据速度优先和压缩比优先算法进行压缩参数自适应调整，在地形分块处理过程中，使用平滑算法解决了边界连续和陡变问题，对于平坦、中等崎岖、高度崎岖的DEM地形数据，本发明分别取得4940倍、1240倍和620倍左右的压缩比。同时，本发明基于CUDA技术对地形解压算法进行了加速，与传统基于CPU的方法相比，极大提高了运算效率，可更好的满足三维地形绘制等实时应用的需要。

Description

一种基于小波自反馈的大规模地形数据压缩和解压方法

技术领域

本发明涉及基于小波自反馈的大规模地形数据压缩和解压方法。

背景技术

三维地形可视化在军事仿真、三维游戏、城市规划、数字地球等领域有着广泛的应用。受计算机内存容量的限制，随着地形规模的增大，需要的存储容量也相应增加。大规模地形由于数据量很大，且应用背景比较复杂，在当前计算机硬件条件下，所需的地形数据往往无法全部载入内存。而且，由于网络带宽的限制，很多需要快速传递地理信息的网络化服务（如防汛、交通等）很难及时完成。战场仿真和三维游戏对地形数据的读取则更为频繁，因而针对大规模地形数据的数据查询和处理效率在这些应用中会受到极大的限制。因此，在内存有限的情况下，如何使载入内存的数据包含更多的地形信息成为虚拟仿真领域一个研究方向，而DEM数据压缩技术则是解决该问题的重要方法之一。

随着小波变换在图像领域的应用日益广泛，利用小波变换对地形数据压缩取了得若干成果，但普遍存在压缩参数选择困难、数据还原效率低下等问题，并且目前的研究大都是在小规模地形上展开，而针对大规模地形压缩技术的研究还相对欠缺。其中，罗永等在数字高程模型数据小波压缩算法中讨论了不同小波基对压缩性能的影响，如提升9-7整数小波、CDF97小波在保持地形形状和起伏特征的前提下对DEM数据最大有4096倍的压缩比。但该方法将所有的地形高程值量均化到0-255范围内，若地形的采样精度很高，该方法会丢失很多的细节。此外，该方法没有说明如何根据地形的变化幅度来动态选择压缩参数。崔宝侠等设计的基于小波的数字高程模型压缩方法，提出了最优标量量化器，可使用该量化器对小波变换后的数据进行自动选择参数量化。但因为其采用简单的游程编码，压缩比并不高，仅具有30倍左右压缩比。李毅等人在2009年提出了基于SPIHT小波的DEM自适应压缩方法，并设计了一种度量地形复杂度的方法，可根据计算出的复杂度对地形的压缩参数进行选择。该方法具有一定的自适应性，但无法精确控制图形的信噪比、压缩比、解压缩的效率等参数。

发明内容

本发明的技术解决问题是：解决了大规模地形数据压缩无法对参数进行自适应选择的问题，并借助基于CUDA的小波变换并行计算解决了压缩算法效率低下的问题，本发明的压缩结果同时具有较大的压缩比和较高的还原效率。

本发明所采用的技术方案为：一种基于小波自反馈的大规模地形数据压缩和解压方法，该方法包括如下步骤：

步骤（1）、大规模地形由于数据量大，很难一次载入内存，对大规模地形数据进行分块预处理，并利用速度优先算法和压缩比优先算法实现了编码过程中对最佳参数的自动选择；

步骤（2）、对地形中的每一个分块，基于小波自反馈的地形数据压缩处理方法来对其进行压缩编码，采用的具有线性相位的双正交9/7小波(CDF97)作为滤波器，具体计算过程中将二维小波变换分解为两次一维离散小波变换来加速计算，即首先将分块后的地形数据和h进行卷积，抽取偶数项，排列作为低频数据；然后将信号和g卷积，抽取奇数项，排列为高频数据；

步骤（3）、由于地形数据的解压缩效率会严重影响到三维地形的实时绘制效率，基于CUDA的小波逆变换方法来加速解压处理：在CUDA架构下，多个并行的thread可以组成一个block，处于同一个block中的thread可共享同一块内存，从而可以进行快速的同步处理，但由于共享内存存储量有限，且只能被一个block中的threads读取，根据CUDA的特点利用专门的数据组织结构来保证算法对中间数据访问的效率；

步骤（4）、通过在分块边界处计算平均高程值的办法来对边界数据进行一致性平滑处理，解决了因为基于小波自反馈的地形数据压缩处理为有损压缩，且分块压缩时由于分块地形在块和块之间缺乏联系，从而解压后会在分块边界处出现数据不一致、过渡不连续的问题。

其中，步骤（1）中所述的速度优先算法（Speed First Algorithm）将期望峰值信噪比、期望像素深度作为比较基准，使得还原后的高程图像满足峰值信噪比高于Epsnr，而每个像素所需要的位数接近Ebpp的要求；压缩比优先算法是在给定期望峰值信噪比Epsnr和小波变换层数Elevel的情况下，自动选择压缩参数bpp，使得还原后的高程图像和原始高程图像的峰值信噪比高于Epsnr，并且具有最高的压缩比。

此外，步骤（3）中所述的基于CUDA的小波逆变换，分为基于CUDA的一维小波逆变换和基于CUDA的二维小波逆变换。

本发明的原理在于：

地形数据压缩算法的基本思想是将地形高程值看作一张二维灰度图，根据相邻地形高程值数据具有连续性的特点，对该灰度图进行小波变换，并用提取低频数据的方法对频域数据进行预处理，以实现高程数据的压缩。压缩算法的处理流程如图1所示：假设原始地形为一个m×n的地形，保存在一个浮点型矩阵O中，矩阵O通过一个量化器HQ将浮点类型的地形数据量化为整数，得到一个量化后的矩阵Q，矩阵Q通过小波变换系统HW将矩阵Q变换为矩阵W。变换需要选择合适的小波基，主要负责将原数据分解为2个部分，即低频平滑部分(LL)和高频细节部分(LH、HL、HH)。其中，子带数据LH、HL、HH分别对应于图像的水平方向、垂直方向及对角线方向的边缘和细节。这个分解过程可重复地用于LL部分的递归分解，重复分解的次数为level。经过变换的矩阵W通过编码系统HC后便可得到压缩后的序列A。

参数bpp(bits per pixel)为每个像素所占的比特数，该参数会直接影响压缩比。对于地形数据来说，每个象素用一个浮点数表示，占用32bit。当bpp=n时，压缩比即为32/n。bpp参数设置的越大，压缩后还原的图像与原始图像的差别就越小。通常采用图像处理中广泛使用的峰值信噪比(PSNR)作为DEM数据的失真测度，来评价还原图像和原始图像的差别。图2给出了一张原始地形灰度图像和bpp=0.026以及bpp=0.013时还原后的图像。从图2中可以看出，对于较大的bpp，图像的PSNR也较大，图像失真较小。地形数据解压算法的流程与压缩算法相反：根据参数bpp，压缩数据A通过反编码系统HC'生成小波矩阵W'，再通过小波逆变换系统HW'生成Q'，最后根据Base值量化成还原数据O'。

本发明与现有技术相比的优点在于：

（1）、针对传统地形压缩算法压缩参数选择困难问题，本发明设计了速度优先算法和压缩比优先算法，实现了编码过程中的对最佳参数的自动选择；

（2）、本发明设计了基于CUDA的小波逆变换方法，提高了解压效率，可更好的满足三维地形绘制等实时应用的需要；

（3）、传统的大规模地形数据压缩方法，在对分块地形数据进行压缩时，因块和块之间缺乏联系会出现数据还原后边界数据不一致、过渡不连续问题，本发明通过设计数据一致性处理方法，解决了这一问题。

附图说明

图1压缩算法系统示意图；

图2原始地形和在不同bpp下的还原效果示意图；其中，图2(a)为原始图像，图2(b)中bpp=0.026，PSNR＝34.1dB，图2中(c)bpp=0.013；

图3边界数据陡变和消除陡变示意图；图3(a)为边界数据陡变问题示意图；图3(b)为消除边界数据陡变后示意图；

图4不同PSNR下还原地形和原始地形比较示意图；图4(a)为原始地形示意图，图4中(b)PSNR=35，bpp=0.1时示意图，图4(c)中PSNR=30，bpp=0.05时示意图，图4(d)中PSNR=20，bpp=0.01时示意图；

图5为PSNR=30时压缩数据还原地形与原地形对比示意图；图5(a)为平坦地形原始数据示意图，图5(b)为平坦地形还原数据示意图；图5(c)为中等崎岖地形原始数据示意图；图5(d)为中等崎岖地形还原数据示意图；图5(e)为高度崎岖地形原始数据示意图，图5(f)为高度崎岖地形还原数据示意图。

具体实施方式

下面结合附图及具体实施方式进一步说明本发明。

本发明中大规模数据的压缩、解压及数据连续性处理算法的实现方法具体如下：

（1）一维小波变化

假设S(0)~S(n-1)为原始信号，为方便卷积计算，将信号扩展到[-4,n+3]上，其中S(-i)=S(i)，S(n-1+i)=S(n-1-i)；(0<i<5)。对于变换后的低频信号a(k)和高频信号d(k)，计算公式如下：

$a\left(k\right)=\underset{i=-4}{\overset{4}{\mathrm{\&Sigma;}}}{(-1)}^{i+1}g\left(i\right)\&times;s(2\&times;k+i)(0\&le;k<n/2)$

$d\left(k\right)=\underset{i=-3}{\overset{3}{\mathrm{\&Sigma;}}}{(-1)}^{i+1}h\left(i\right)\&times;s(2\&times;k+1+i)(0\&le;k<n/2)$

将a(k)，d(k)排列成a(0)~a(n/2-1)；d(0)~d(n/2-1)的形式即实现了一维小波的正变换。

小波变换的逆变换和正变换类似，但过程相反：先将a(0)~a(n/2)和d(0)~d(n/2)扩展到a(-2)~a(n/2+1)和d(-2)~d(n/2+1)上，其中a(-i)=a(i)，a(n/2-1+i)=a(n/2-i)，d(-i)=d(i-1)，d(n/2-1-I)=d(n/2-1-i)；（0<i<3）；然后，分别将其与h，g进行卷积即可得到原始信号S(0)~S(n-1)。

（2）基于CUDA的一维小波逆变换

假设信号a(0)~a(n/2-1)；d(0)~d(n/2-1)保存在全局内存中，每一个block包含n/2个线程，为了不直接读取全局存储器，定义两个长度为(n/2+4)的共享的内存sa，sd。这两个共享内存的寻址范围为[-2,n/2+1]。此外还定义一个长度为n的共享存储器ss，用于保存结果。算法并不直接将计算结果回写到全局内存中，而是先写到共享的内存ss中，具体流程如下：

1）当前线程的索引为tID，用两条复制语句sa(tID)=a(tID)，sd(tID)=d(tID)进行复制，该访存操作满足全局内存的访问条件，可并行完成；

2）同步同一block中的所有线程，并用前两个线程将sa，sd扩展到[-2,n/2+1]中；

3）同步block中所有线程，运用卷积公式计算ss(tID*2)和ss(tID*2+1)的值，其中ss(tID*2)=h(0)×sa(tID)+h(2)×((sa(tID+1)+sa(tID-1))+g(3)×((sd(tID+1)+sd(tID-2))+g(1)×((sd(tID)+sd(tID-1))。该公式对于sa，sd的访问满足共享内存的访问条件，可并行完成，而对ss的写入需要互斥进行；

4）同步block中所有线程，将ss中的数据回写到全局内存中，该访存操作满足全局内存的快速访问条件，可以并行完成。

（3）基于CUDA的二维小波逆变换

基于CUDA的二维小波逆变换算法为：假设待解压数据W是一个N×N的矩阵，首先对W中每一列进行一维小波逆变换，变换后的矩阵为P，然后对矩阵P中每一行进行一维小波逆变换，变换后的矩阵为Q。

由于CUDA中对列的访问会带来大量的内存冲突，因此本发明采用了求矩阵转置的办法把列变换变为行变换。实现流程如下：

1）将矩阵W从内存复制到全局内存中；

2）对矩阵W进行求转置运算，得到矩阵IW；

3）使用N个grids，每个grid里N个threads，采用一维小波变换对IW中每一行进行变换，并且每一行分别需要一个grid来进行单独处理，得到变换后的矩阵IP；

4）对矩阵IP进行求转置运算，得到矩阵P；

5）使用N个grids，每个grid包含N个threads，对P中每一行进行变换，得到最终结果Q；

6）将结果从GPU全局内存中复制到内存中。

（4）数据边界一致问题处理

设初始的大规模地形被分成了K行，L列，这里用B(I，j)表示一个块，在地形数据被压缩和解压后，数据是一个N×N的矩阵O，用O(I，j)(0≤I，j<N)表示地形块中第i行、j列点的高度值。将矩阵O扩展成(N+1)×(N+1)的矩阵R，即R(I，j)=O(I，j)，R(N,i)=O(N-1,i)，R(i，N)=O(I,N-1)，R(N,N)=O(N-1,N-1)；(0≤I，j<N)。

在原始地形中的一个点可能同时属于多个不同的块，如第B(i，j)块中的R(i,N)和第B(i，j+1)块中的R(i,0)描述的是同一个点。对每个块进行压缩解压处理，会使得不同块中同一个边界点的高度出现很小的不一致，这种不一致在地形绘制时出现裂缝现象。本发明采取计算平均值的方法来处理边界数据，处理方法如下：

1）如果一个边界点a，仅仅属于一个块，则a的高度处理后和在该块中对应点的高度一致；

2）若数据a同时属于两个块，不论a是边界数据或是顶点上的数据，在进行数据解压缩时，a的高度取a在这两个块中所对应数据的平均值；

3）若数据a为顶点数据，且同时属于4个块，则进行数据解压时a取这4个块中对应数据的平均值。

通过上述处理后，块和块之间边界数据不一致的问题可以得到解决，但如图3(a)的示意图所示，该操作会导致边界数据波动过大，出现类似鸿沟的地形，我们进一步采用数据平滑处理方法来进行解决。

（5）一维数据平滑处理

平滑的定义：设x(i)和x(i+1)为任意连续的两个数据(n-1<i<m)，若对于任意的i，abs(x(i)-x(i+1))<δ则称x(n)到x(m)为平滑的，其中δ为门限值。

对于一维数据x(0)~x(N)，若x(1)~x(N-1)是平滑的，而x(0)~x(1)和x(N-1)~x(N)不平滑，则通过如下算法将x(0)~x(N)变为平滑的一维数据y(0)~y(N)：

1）两端数据不变即y(0)=x(0)，y(N)=x(N)；

2）将x(0)~x(1)和x(N-1)~x(N)之间的高度差分等分为N/2份，每一份为t1=(x(1)-x(0))×2/N和t2=(x(N-1)-x(N))×2/N；

3）令x(1)=x(0)，x(N-1)=x(N)；

4）将高度差均匀的分配到各个连续的数据上，即y(i)=x(i)-t1×(n/2-i)i<n/2；y(i)=x(i)-t2×(i-n/2)i>n/2。

通过以上步骤数据x(0)~x(N)变为了平滑的一维数据y(0)~y(N)。

（6）二维平滑处理

对于一个(N+1)×(N+1)的矩阵R，二维地形的平滑处理过程如下：

1）将处于地形四个边界的数据运用一维数据平滑处理算法进行处理，使得边界平滑；

2）将地形数据的每一行运用N次一维数据平滑处理算法进行处理，使得每一行平滑；

3）将地形数据的每一列运用N次一维数据平滑处理算法进行处理，使得每一列平滑。

经过处理后的地形如图3(b)的示意图所示，消除了鸿沟，具有良好的平滑连续性，且保持了原地形的地貌特征。

本发明在具有1.6GHz CPU，GeForce 8600GS GPU，2GB RAM的PC机上进行了实验，实验数据为中国内蒙某地60km×70km的真实地形数据，数据采样精度为2.4米，高程值的误差Δ=0.5m。

实验中将地形数据分成112×96共10752块，每一块包含256×256个高程值。我们主要进行了两个方面的实验：一是随机选取一块地形，测试该地形对于不同峰值信噪比的压缩倍数，以及解压后的地形与原始地形相似的程度；二是分别对平坦地形、中等崎岖地形和高度崎岖地形进行压缩和解压测试，测试在失真度较小的情况下所能达到的最大压缩倍数，实验结果请见表1和表2。

表1压缩比随峰值信噪比的影响

表2PSNR=30时不同崎岖程度地形最高压缩比

表1列出了地形压缩比随峰值信噪比的变化情况，由表1可以看出，压缩比随着峰值信噪比的升高而降低。如：图4给出了PSNR=35、30、20时原始地形和还原地形的效果对比示意图，在PSNR=30时，还原地形的失真度较低。所以可以在PSNR=30下来测试最高压缩比。表2为对三块随机选取的平坦地形、中等崎岖地形和高度崎岖地形下最高压缩比的测试数据。从表2可以看出，在PSNR=30的情况下，压缩比可分别达到4940、1240、624；而且三块地形所能达到的最大压缩比随着地形崎岖程度的增加而减小，此外，图5给出了PSNR=30时原地形和还原地形的效果对比示意图。

本发明未详细阐述的技术内容属于本领域技术人员的公知技术。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (3)

1.一种基于小波自反馈的大规模地形数据压缩和解压方法，其特征在于：该方法包括如下步骤：

步骤(1)、大规模地形由于数据量大，很难一次载入内存，对大规模地形数据进行分块预处理，并针对速度优先算法和压缩比优先算法实现了编码过程中对最佳参数的自动选择；

步骤(2)、对地形中的每一个分块，基于小波自反馈的地形数据压缩处理方法来对其进行压缩编码，采用的具有线性相位的双正交9/7小波(CDF97)作为滤波器，具体计算过程中将二维小波变换分解为两次一维离散小波变换来加速计算，即首先将分块后的地形数据和h进行卷积，抽取偶数项，排列作为低频数据；然后将信号和g卷积，抽取奇数项，排列为高频数据；

步骤(3)、由于地形数据的解压缩效率会严重影响到三维地形的实时绘制效率，基于CUDA的小波逆变换方法来加速解压处理：在CUDA架构下，多个并行的线程(thread)可以组成一个线程块(block)，处于同一个线程块(block)中的线程可共享同一块内存，从而可以进行快速的同步处理，但由于共享内存存储量有限，且只能被一个线程块(block)中的多个线程(thread)读取，根据CUDA的特点提出专门的数据组织结构来保证算法对中间数据访问的效率；

步骤(4)、通过在分块边界处计算平均高程值的办法来对边界数据进行一致性平滑处理，解决了因为基于小波自反馈的地形数据压缩处理为有损压缩，且分块压缩时由于分块地形在块和块之间缺乏联系，从而解压后会在分块边界处出现数据不一致、过渡不连续的问题。

2.根据权利要求1所述的基于小波自反馈的大规模地形数据压缩和解压方法，其特征在于：步骤(1)中所述的速度优先算法(Speed First Algorithm)将期望峰值信噪比、期望像素深度作为比较基准，使得还原后的高程图像满足峰值信噪比高于期望峰值信噪比Epsnr，而每个像素所需要的位数接近每个像素所占的比特数bpp的要求；压缩比优先算法是在给定期望峰值信噪比Epsnr和小波变换层数level的情况下，自动选择压缩参数bpp，使得还原后的高程图像和原始高程图像的峰值信噪比高于Epsnr，并且具有最高的压缩比。

3.根据权利要求1所述的基于小波自反馈的大规模地形数据压缩和解压方法，其特征在于：步骤(3)中所述的基于CUDA的小波逆变换，分为基于CUDA的一维小波逆变换和基于CUDA的二维小波逆变换。