基于协方差消除多重共线性的一种方法

孟金玉 袁永生

摘 要：针对传统三维地形生成算法在生成大规模地形数据时用时较多的问题，提出了一种结合双线性插值和Perlin Noise的地形生成算法。该方法首先将初始地形数据进行扩展，之后利用双线性插值和Perlin Noise进行地形细节生成。该方法能够有效降低生成大规模地形数据时的时间，并且该算法还可利用低分率地形数据生成高分辨率地形数据。

关键词：三维地形; Diamond-Square; 双线性插值; Perlin Noise

中图分类号：TP391      文献标识码：A

三维地形生成技术一直是一个热门的研究方向，國内外研究人员对此进行了大量的研究，并提出许多生成方法。目前三维地形的生成有多种形式，其中研究最为广泛的有两种，第一种是基于程序的地形生成方法，一般采用基于分型几何的方法生成;另一种是基于数据的地形的生成方法，一般为基于数字高程模型的生成。

分形几何算法主要有泊松阶跃法、中点位移法、逐次随机增加法、傅里叶滤波法、带限噪音累计法等五种方法[1]。其中，中点位移法是最基础和应用最广泛的方法，较为著名的有Miller提出的“Diamond-Square”算法[2]，该算法在生成大规模高程数据时，存在时间较长的问题。

提出一种新的地形生成算法——Bilinear-Perlin算法。该算法将双线性插值与Perlin Noise算法有效融合，通过在低分辨率高程数据基础上进行插值获得高分辨率数据。利用该算法能够有效降低“Diamond-Square”算法在生成高分辨率大规模地形数据时的生成时间，并能够在低分辨率的数字高程模型下通过插值获得高分辨率地形。

1 Diamond-Square算法

Diamond-Square（以下简称DS）算法主要分为正方形和菱形两个步骤。迭代这两个步骤，不断计算获得高程数据，细化地形内部网格。

算法如图1所示，现在作如下定义：Size表示欲生成高程数据的数据规模边长，例将生成的数据规模为512*512，则Size为512。RectSize表示当前步长，此变量在一次正方形步骤后，缩小为原来的1/2。已知A、B、C、D四点数据，分别记为Ha、Hb、Hc和Hd。

由DS算法效果图能够看到，DS算法生成的效果较为真实，但是由于算法特性，在生成小规模数据时速度较快，生成大规模数据时，速度较慢。

3 Bilinear-Perlin算法

3.1 算法思想

Bilinear-Perlin（以下简称BP）算法的核心思想为插值计算，即通过对小规模高程数据进行扩展，扩展到大规模数据后，产生的空白区域通过插值计算获得。空白区域数据首先采用双线性插值获得，但由于双线性插值算法得到的数据较为光滑，与真实地形不符，通过Perlin Noise算法添加随机值，增加地形表面粗糙度，以此得到较为真实的高程数据。最后，利用平滑滤波函数对所有数据进行一次平滑滤波，最终得到插值计算后的大规模数据。BP算法流程图如图5所示。

3.2 双线性插值

双线性插值是有两个变量的插值函数的线性插值扩展，其核心思想是在两个方向分别进行一次线性插值[5]。

三维地形高程数据通常存放于二维数组中，数组坐标表示地形点的屏幕坐标，数组的值表示当前点的高程值。如图6所示，在已知点（0，0）、（0，1）、（1，0）（1，1）四点值，求点（x，y）的值，其中0≤x≤1，0≤y≤1。此时需要利用双线性插值获得点（x，y）的值。

对于计算任意给定的（x，y）坐标上的精确高度，首先计算x和y值，以xRelative和yRelative表示。代码如下：

int xLow= xCoord;

int xHig=xLower+ 1;

float xRelative = （xCoord-xLow） / （xHig- xLow）;

int yLow= yCoord;

int yHig=yLower+ 1;

float yRelative= （yCoord-yLow） / （yHig- yLow）;

在获得xRelative和yRelative后，需要获得其所对应的具体高程值。三维地形的绘制采用三角形的方式进行绘制，但是可以有两种方式进行绘制，如图7所示。两种地形绘制方式的不同，会影响到P点的高度。

在地形绘制中，通常采用右侧的绘制方式，因为右侧的绘制方式能够更容易的确定点P（x，y）在哪个三角形上。在右图中，如果xRelative + yRelative等于1的话，P（x，y）在对角线上;如果和小于1，P（x，y）位于左下三角形内;否则，P（x，y）在右上三角形内。

在获得了对应点的坐标，以及点位于哪个三角形中，通过四个顶点的高度利用双线性插值的方法获得精确高度。

定义：boolpointAboveLowerTriangle。

如果点在左下三角形中，pointAboveLowerTriangle为true，否则为false。

其中点在左下三角形中，计算方法如下：

finalGht = ghtLxLy;

finalGht += yRelative*（ghtLxHy-ghtLxLy）;

finalGht += xRelative*（ghtHxLy-ghtLxLy）;

点在右上三角形中，计算方法如下：

finalGht = ghtHxHy;

finalGht += 1.0f-yRelative * ghtHxLy-ghtHxHy;

finalGht += 1.0f-xRelative * ghtLxHy-ghtHxHy;

其中，LxHy表示在四个点中，x值较低与y值较高的坐标，即左上角坐标，其他坐标同理。根据以上计算方法，我们就能够获得所有点的插值数据，但是这些点的数据由于采用线性计算，插值产生的地形过于平滑，与真实的地形不符，因此我们需要进行添加随机高度，使地形产生粗糙表面。

3.3 Perlin Noise

1985年，Ken Perlin首次发表了Perlin Noise[6]算法。Perlin Noise算法通过混合不同频率不同振幅的噪声函数，使得相邻数据差异较小，使整体数据符合自然环境的随机状态。目前，人们通常使用Perlin Noise算法生成木头、石头、云状和火焰等纹理。

Perlin Noise算法由两部分组成，噪声函数和插值函数。噪声函数用来随机产生整点数据，插值函数用于平滑整点数据。

3.3.1 噪声函数

噪声函数的本质为一个随机生成函数，其是离散的。通过输入整数参数，返回一个随机数。如图8所示，输入x值，对应产生一个0到1之间的随机数，但是当输入的x值相等时，输出的随机数相等，此特性为Perlin Noise算法的最主要特性。由于其是离散的点组成的，因此我们只需要构建离散点数据即可，并且噪声函数的频率和振幅都能自行控制。

3.3.2 插值函数

噪声函数产生的值是离散的，因此需要通过插值函数，填补两个离散点之间的部分，生成一个非整数参数的连续函数。目前常用的插值函数共有三种：

（1）线性插值函数：

A×1-x + B×x（7）

参数：A，B为需要插值的两个离散点终点。x取值为0到1，表示插值在两点间的比例。x=0表示在A点，x=1表示在B点。

（2）余弦插值函数：

fx=1-cos（xπ）×0.5

return A×1-fx +B×fx（8）

参数：同线性插值函数。但是此种插值函数更为平滑。

（3）立方插值函数：

P=（v3-v2）-（v0-v1）

Q=（v0-v1）-P

R=v2-v0

S=v1

return Px3 + Qx2 + Rx + S（9）

参数：立方插值函数需要五个参数v0、v1、v2、v3、x，v1，v2，x同上面参数，v0，v3为前一点和后一点。此种方法最为平滑，但是效率较低。

3.3.3 Perlin Noise函数

在介绍Perlin Noise函数之前，首先需要明确两个概念，振幅（Amplitude）和频率（Frequency）。如图9所示，振幅表示最大值和最小值的差值，频率表示为1/波长。

如图10示不同频率和振幅的噪声函数。可以看到，通过改变振幅和频率，能够得到不同样式的噪声函数。

Perlin Noise的核心就是将不同频率和振幅的噪声函数进行叠加，如图11所示，叠加了图10所示的噪声函数。大振幅低频率的噪声函数决定了图形的形状，例如山脉的轮廓;小振幅高频率的噪声函数决定了图形的细节，例如地形的侵蚀状。

以上分别介绍了双线性插值和Perlin Noise算法，每次利用双线性插值算法获得插值数据后，添加相应的Perlin Noise函数点。这样插值获得的数据不仅能够和原有地形走势相同，并且还能够与真实的地形一样，产生地形的粗糙感，使其更真实。

4 DS算法改进

前文已经介绍，由于DS算法的原因，其在生成大规模地形数据时，生成时间会非常大。下面利用BP算法对DS算法进行改进，首先采用DS算法生成小规模数据，然后利用BP对小规模数据进行插值扩展为大规模数据。此方法能够在保证较高地形细节的基础上，生成速度得到明显提高。

将改进的DS算法称为DS-P算法。DS-P算法的算法流程图如图12所示。

5 算法比较

5.1 DS-P算法与DS算法比较

由图[4]和图[13]对比可知，DS-P和DS算法在生成地形的趋势和细节上基本是一致的，DS-P算法生成的地形满足基于程序的地形生成算法的要求。

下面将对两种算法生成时间进行对比。实验环境如下所示。

1.硬件环境

CPU： Intel Core i3-2120  主频3.30GHz

内存：4GB

显卡：Intel Graphics Family

2.软件环境

操作系统：32位Windows7 SP1

编译环境：VS2010

编译语言：C#

实验采用多次平均方法，在同样环境下，各种算法分别进行生成地形数据，各生成50次，然后平均处理，得出生成时间，其中DS-P算法的初始规模数据位64×64。时间统计如表1所示。

图15为两种算法所用时间的曲线图。由图可以看到随着地形数据规模的增加，DS算法的地形生成所需时间会陡然加剧，而改进的DS-P算法耗时会远远小于DS算法。

5.2 BP算法生成高分辨率地形

BP算法不仅能够有效降低DS算法在生成大规模地形数据时的生成时间，并能够在低分辨率的数字高程模型下通过插值获得高分辨率地形。

现将三种高程数据进行对比，第一组：原始数据分辨率30m×30m;第二组：在第一组的基础上进行间隔采样，每三个点采样一个，分辨率降低为90m×90m;第三组：利用BP算法对第二组数据进行插值，使分辨率重新达到30m×30m 。渲染显示如图16所示。

由上图可知，中间图形与左图对比，地形出現明显折痕，细节度几乎没有，对原有数据的还原效果极差。右侧图形与中间图形对比，细节度明显提高，与左侧图形对比，差异较小，如图17所示。

6 结 论

DS算法在生成大规模地形数据时，会存在生成时间过长的问题。本文提出BP算法，并根据此算法提出了改进的DS-P算法。通过实验表明，在不改变原有地形细节时，DS-P算法在生成相同地形数据规模时，所用时间大大减少。并且，BP算法能够应用到将现有的低分辨率的数字高程模型下通过插值获得高分辨率地形，细节还原较高。

参考文献

[1] MUSGRAVE F K ， KOLB C E ， MACE R S . The synthesis and rendering of eroded fractal terrains[J]. Computer Graphics， 1989， 23（3）：41-50.

[2] MILLE R， GAVIN S P. The definition and rendering of terrain maps[J]. Computer Graphics， 1986，20（4）：39-48.

[3] 朱庆.分形理论及其在数字地形分析和地面仿真中的应用[D]. 北京：北方交通大学，1995.

[4] 雷磊. 三维地形生成及可视化技术研究[D].哈尔滨：哈尔滨工程大学，2005.

[5] 叶金印， 邱旭敏， 黄勇，等. 气象遥感图像及格点场重采样插值方法[J]. 计算机工程与应用， 2013， 49（18）：237-241.

[6] PERLIN K. An image synthesizer [J].Computer Graphics，1985，19（3）：287-296.