CN105336003B

CN105336003B - 结合gpu技术实时流畅绘制出三维地形模型的方法

Info

Publication number: CN105336003B
Application number: CN201510625734.3A
Authority: CN
Inventors: 全吉成; 王玉锋; 杨明权; 白新伟; 孙文邦; 王平; 赵柏宇; 候宇青阳; 魏湧明
Original assignee: AIR FORCE AVIATION UNIVERSITY OF CHINESE PEOPLE'S LIBERATION ARMY
Current assignee: AIR FORCE AVIATION UNIVERSITY OF CHINESE PEOPLE'S LIBERATION ARMY
Priority date: 2015-09-28
Filing date: 2015-09-28
Publication date: 2018-05-25
Anticipated expiration: 2035-09-28
Also published as: CN105336003A

Abstract

一种结合GPU技术实时流畅绘制出三维地形模型的方法，属于影像处理技术领域。本发明的目的是以全球数字高程模型为数据源，基于当前热门的可编程GPU技术，实现了多次绘制时的缓存复用，有效地减小了计算空间负载的结合GPU技术实时流畅绘制出三维地形模型的方法。本发明的步骤是：多分辨率金字塔模型构建、影像噪点的消除、对影像进行滤波、对平面投影后的地球按照等经纬度进行分块，同时按照自顶向下的方式，构建不同层次级别的金字塔层。本发明基于可编程GPU技术实现了地形渲染的加速和增强，即通过使用着色器语言控制图形绘制流水线的各个阶段，将高程数据的顶点信息和索引信息分别生成两张纹理存入显存，供整个地形绘制调用；利用曲面细分和分形技术在几何阶段对顶点插值和偏移，生成过程化细节，弥补分辨率不足时地形网格棱角化的现象。

Description

结合GPU技术实时流畅绘制出三维地形模型的方法

技术领域

本发明属于影像处理技术领域。

背景技术

随着航空航天科技和遥感技术的日益成熟，环绕在近地轨道的侦察卫星以及各种飞行器上所携带传感器的空间分辨率在不断提高，数字摄影测量技术也越来越精细准确。多种遥感数据所汇聚而成的数据集是地理空间情报的重要信息来源，而其中的遥感影像和数字高程则是整个空间数据框架的基础和根基。在对这种基础地理数据的应用需求急速增长的前提下，卫星遥感影像和数字高程的数据量也跟着日益庞大起来，如何对之进行有效管理和利用，成为了遥感数据处理和地理信息系统领域里迫切需要解决的问题，特别是发布管理和组织调度海量的遥感信息，在为国土规划、测绘遥感、军事交通、水利电力、农林环保等各种应用需求提供地理信息服务上具有极其重要的现实意义。

现代军事技术的发展，正在不断把战场空间向着多维度、多要素的方向推进。面对复杂交错、瞬息万变的战场环境和态势，如何协助指挥员全面地部署作战行动、准确地掌握战场信息就成了军事情报的重要目标和内容。与此同时，我们的作战空间在不断扩展，而地理空间情报的发展却略显滞后，早期由于传感器技术和计算机图形技术的发展限制，数据的获取和可视化能力难以达到实时准确的要求，二维平面的处理方式可以简化空间信息理解和表达的过程，但却是以牺牲空间信息的真实性和完整性，尤其是高程信息和三维空间的拓扑关系为代价的。然而当遥感科学技术和计算机图形学的日新月异，使获取高分辨率的高程和影像，高效逼真的三维显示不再是难以实现的技术挑战，学科领域的发展也都促使军事情报必须弥补其三维表达能力不足的技术缺陷。

在军事领域，战争形态的发展要求指挥员对战场环境的掌握不能再局限于某一固定区域或某种特定形式的情报信息上，想要获取现代信息化条件下局部战争的胜利，就必须及时灵活、全方位地掌控战场，而能否实时准确清晰全面地展现战场环境的诸多要素，为指战员进行辅助决策，是军事情报面临的严峻挑战。

发明内容

本发明的目的是以全球数字高程模型为数据源，基于当前热门的可编程GPU技术，实现了多次绘制时的缓存复用，有效地减小了计算空间负载的结合GPU技术实时流畅绘制出三维地形模型的方法。

本发明的步骤是：

三维地形模型的多分辨率金字塔模型构建：

①地形数据的预处理流程：对原始影像进行滤波得到无噪声的影像，并构建多分辨率金字塔；

Ⅰ、影像噪点的消除：

噪声模型符合加性噪声的特征：

其中，I(x,y)为位于影像(x，y)处的灰度值，O(x，y)为未受噪声污染的灰度值，n(x，y)是独立于O(x，y)的随机噪声值，p为受到噪声污染的概率；

对影像进行滤波的步骤流程如下：

(ⅰ)、计算出整幅图像的均匀度image_HSV，然后对整幅图像的像素进行遍历；

均匀度是衡量影像中噪声信号对比度的参数，影像I_M×N和窗口W_r×r的均匀度分别定义如下：

假设图像I大小为M×N，影像中任一点(x，y)的灰度值为f(x，y)，整幅图像的平均灰度值average(I_M×N)的计算公式为：

滑动窗口的大小为r×r，若窗口中心像素的坐标为(i，j)，灰度值为f(i，j)，则该窗口内的像素平均值average(W_r×r)的计算公式为：

(ⅱ)、统计当前窗口内的极大值点像素值iMax，及其个数nMax，窗口内的像素个数sumPixel，以及窗口的块均匀度wnd_HSV；

(ⅲ)、判断当前窗口中心像素点(x，y)的灰度值I(x，y)是否等于窗口内像素极大值iMax，或者窗口均匀度wnd_HSV是否大于影像均匀度image_HSV；

(ⅳ)、对判定为噪声点的四个方向进行拉格朗日插值；拉格朗日的插值方法为：

其中，y_i表示x_i出的函数值；此处，插值算法是针对一维数据进行处理，当拓展到二维影像上时，取插值像素点(x，y)水平、垂直、两个对角线方向共四个方向上的临近像素点进行插值，并取其平均值作为最后的插值结果；

②通过上述清除噪点后的数据构建多分辨率金字塔：

金字塔结构是一种多分辨率层次模型，采用倍率方法构建，即对同一地区，采用一系列的网格表示，上下相邻网格的采样间隔相差相同的比例系数m，这样从金字塔自底向上，分辨率逐层降低，网格尺寸也相应减小；

设原始网格的尺寸为S₀，分辨率为R₀，所构建的金字塔共有l层，金字塔顶端为第0层，那么第k层的网格尺寸s和分辨率r可以表示为：

对平面投影后的地球按照等经纬度进行分块，同时按照自顶向下的方式，构建不同层次级别的金字塔层，构建的方案如下：

(Ⅰ)地球的纬度跨度为-90°到90°，经度跨度为-180°到180°，在投影平面上的宽高之比恰为2:1，因此四叉树在第0级时即分为8个根节点，横向纵向跨度均90°，以满足等经纬度划分的要求，除此范围以外的坐标值均视为无效；

(Ⅱ)每层LOD的分辨率之比为2，即子块的分辨率是其父块的二倍，相同的地理空间，子块的像素尺寸是上一级的2倍，而在相同的屏幕分辨率下，父块覆盖的地理空间面积是子块的4倍；

(Ⅲ)每层的块序列采用简单的行列编码，下一层与上一层之间的块数之比、每层横向与纵向的块数之比均为2:1，按照计算机图像坐标系的习惯，起点定义在左上角，行编号由上至下，列编号从左到右；第0层的分块数是4×2，第k-1层的分块数是2^k+1×2^k；

(Ⅳ)每个分块中包含了对应地理位置的DEM数据，其分辨率尺寸必须是奇数，即(2n+I)×(2n+I)；

(Ⅴ)金字塔的底层为最高分辨率，直接从原始数据中获得，其所在的LOD层级由其自身分辨率与所属范围决定；每个地形块对应的纹理数据和法向数据分辨率必须要高于DEM的分辨率，以保证插值平滑，避免渲染失真；假设某一层DEM分辨率为r，这里把它们的分辨率定义为r′=n×(r-l)，n＞l，n∈N；得到了全球规模的金字塔模型。

基于GPU的流水线绘制增强技术，上述全球规模的金字塔模型的数据进行绘制：

利用着色器语言将高程的顶点信息和索引信息分开存储，并作为纹理数据一次性传入显存；

a、节点的数据结构与组织：

采用Chunked LOD算法构建了块状四叉树结构，即不再以单个顶点或三角形为处理单元，而是以包括纹理贴图在内的网格模型，将三角形打包成条带一次性完成渲染，充分利用GPU的批处理能力；

b、缓存的复用：采用Geo-Clipmap所使用的方法，同一层中的顶点坐标(x，y)均进行同样的变换；

c、金字塔模型的数据实时绘制：分为两类，一种是静态节点，另一种是动态节点；

Ⅰ、静态节点数据渲染：即在预处理过程中通过采样获取同一区域不同分辨率节点，在渲染时，通过将对应分辨率的节点调入显存中直接进行渲染；

①环状取址的位运算：

使用一个固定大小的网格缓存获取四叉树中任意节点的高程数据，通过设定偏移地址和地理空间跨度，渲染整个地形空间；在加入视椎裁剪之前，取视点垂直向下，视点的可视范围是一个n×n固定大小的像素网格，根据视点所在的金字塔层级l，得到该层中每个像素所代表的地理空间分辨率r_l＝2^-l-9×π×R，其中R表示地球半径，用(gx_diff，gy_diff)表示视点投影在地形平面上的位移，那么视点在平面上移动的像素距离表示为和视场所覆盖的像素范围由(x_min，y_min)和(x_max，y_max)表示；

移出视场的节点被第一次进入视场的节点代替，其余节点不动，这样就完成了增量更新，对于视场范围内任意一点(x，y)，其在节点数组中存储的位置为(x mod n，y modn)，其中“mod”表示求模运算；

②误差度量的确定：

采用Chunked LOD的误差度量计算方式，即从叶子节点开始，沿着四叉树向上逐层计算节点的误差，一个节点的误差度量ε_c为其四个子节点的ε_a的最大值加上其本身的误差度量ε_■；ε_c的计算法如下：

这样，屏幕空间误差大于阈值时，选择继续向下寻找误差更小的节点，直至满足要求或者到达叶子节点；

③可见性剔除：

视点较近时的剔除方法：

当视点距离地形较近时，采用包围盒与视椎体进行相交测试，即通过包围盒的特征点与视椎体的6个平面依次求交，以确定这些特征点是否属于平面正向的那半个空间或者与平面相交；视椎的6个平面由方程f(x，y，z)＝ax+by+cz+dw＝0决定，其中向量(a，b，c)表示该平面的法向量，数值dw表示平面与原点的距离；包围盒的9个特征点由8个角点和一个中心点组成，通过比较包围盒的特征点即可判断地形网格的可见性，取size等于包围盒的斜对角线长的一半；

视点较远时的剔除方法：

当视点远离地形时，利用节点与视点的距离系数进行判断，定义距离系数大小为节点包围盒所代表的实际距离，在投影平面长度L与包围盒距投影平面上的距离D之比，即a＝L/D；

取视点坐标v(x₀，y₀，z₀)，包围盒中心坐标为c(x₁，y₁，z₁)，则视点到地面的方向向量为d(m，n，p)＝v-c；

设包围盒的8个顶点坐标为p_i(x_i，y_i，z_i)，过点P_i与d平行的直线方程为：

f(x，y，z)＝m(x-x_i)+n(y-y_i)+p(z-z_i)＝0 (3.9)

视椎体的近剪切面作为投影平面，与方向向量d垂直，设其平面方程为

mx+ny+pz＝0 (3.10)

联立上述两个方程即可求得p_i在上述平面上的投影点坐标p_i′为

计算出所有8个顶点的p_i′后，取8个坐标点中的极大极小值min X，max X，min Y，max Y，min Z，maxZ；令

距离系数a的值为：

判定时，在视点高度大于设定阈值时，若节点包围盒的a值小于1，即在投影平面上大小不足一个像素时，即可认为不可见，不再渲染，否则继续渲染该节点；

④缝隙消除与几何过渡：

在进行节点筛选前，首先计算每个LOD层级所覆盖的地理范围，以便渲染时在视窗平面上产生颗粒度近似一致的三角形数量；

设每层LOD的起始距离为start_i，i表示当前LOD所在的层级，其所覆盖的距离范围为range，而视点到当前LOD任意一个顶点的距离为dist，则当前的过渡区域[mStart,mEnd]以及该顶点的过渡系数k为：

其中，clamp(x，a，b)为钳位函数，即

获得过渡系数后，定义顶点所在网格尺寸为n×n，所在的四叉树深度为l，顶点在网格内归一化的坐标为(x，y)，则网格尺寸相对于世界坐标下的实际距离放缩值scale以及过渡后的坐标(x′，y′)由如下公式求得：

其中，frac(x)表示取x的小数部分；

最终，顶点坐标会从(x，y)变换到(x/2，y/2)的位置，它对应的高程通过双线性插值获得；

Ⅱ、动态节点：是在四叉树结构延伸至叶子节点时，仍不能满足视点阈值的需求，由GPU实时计算生成并加入四叉树中的节点；

动态节点的生成：

①地形模型的细分：细分是在四边形网格的原始顶点上根据细分掩膜插入面点和边点，并与原始顶点连接成新的四边形，生成的非边界顶点度数为4，边界顶点度数为3，都是正则点；

正则点细分过程：对四边形V₅V₆V₇V₈进行细分时，会在其四条边上形成四个边点E₀E₁E₂E₃，同时还有一个面点F，使原来的面片分裂成四个子面片；E₀，E₁，E₂，E₃的值如下公式：

面点F的值由公式4.2得出：

其中，

②随机分形：用变量(x，y)表示曲面上点的坐标，X(x，y)表示该点在曲面上的高；则FBM曲面定义为在某概率空间上，指数为H(0＜H＜1)的一个随机过程

随机变量X：R→R²满足：

(1)以概率1成立X(0，0)＝0，即过程自原点开始；且X(x，y)是关于(x，y)的连续函数；

(2)对任意变量(x，y)，(h，k)∈R²，其二维增量X(x+h，y+k)服从期望为0，方差为(h²-k²)^H的正态分布，其概率满足

采用两种分形算法来生成不同样式的噪声图，分别是Diamond-Square算法和Perlin噪声：

Diamond-Square算法：在迭代过程中每一次细分都会执行两步操作，Diamond步和Square步；Diamond步用于生成正方形的中心，取四个顶点计算平均值，再加上一个随机量，作为中心点的值；Diamond步生成了四个三棱锥；Square步则是取棱锥的四个顶点，用其平均值加上随机变量作为原正方形的四条边中点，每一次细分的结果一致；

Perlin噪声：对于一幅噪声图像，首先给整数坐标点赋予随机方向的法向量与灰度值，除此以外任意一点的像素值由该点所在最邻近的单元格内四个顶点的法向量，和连接该点与四个顶点的方向向量来决定，设为图中所求点的向量，n为噪声的叠加次数，一般取n∈[3，12]；设噪声函数波长为λ，则频率ω＝λ^-0.5i，该点的灰度值为：

其中，offset为偏移量取-0.1，是一个噪声函数；

采取预处理的方法生成多种类型的噪声纹理，避免实时计算增加GPU的额外开支，在光栅化阶段中，像素着色器读取一张包含噪声的纹理图与高程图进行融合，同时噪声图进行不同比例地放缩与叠加，扩展出无限大的噪声纹理，以此消除绘制出的周期性和模式化效果。

本发明的“数字地球”和“信息格网”的提出，为解决以上问题提供了新的方法和手段，即对空间数据进行高效地组织管理，并以数字化、可视化的形式表达，而本文所研究的地形信息的可视化，正是将战场环境可视化首要而基本的内容：应用上，任何作战都是在一定的战场地理空间中进行的；技术上，地理空间环境的可视化是其他要素可视化叠加显示的基础和载体，对研究军事情报的三维表达有着重要的意义。其优点是：

1.分析了地形高程的建模方法和基本特点，探讨了该领域的研究现状和热点问题，确立了本文的技术路线，分析解决了高程数据的预处理问题如平面投影坐标系的变换和高程图中噪点的消除等，为组织和管理高程数据奠定了基础。

2.提出了使用金字塔模型构建全球多分辨率虚拟地形环境，并结合多分辨率LOD模型技术，设计并实现了全球范围的金字塔模型，包括金字塔的构建规则和分层分块方案，采用四叉树结构对数据进行索引，利用位运算和队列结构实现数据节点快速定位和更新。

3.分析比较了几种经典的地形LOD算法，结合Chunked LOD算法的四叉树结构和Geometry Clipmap算法的环状取址，实现全球的地形绘制和渲染模型的增量更新；提出了基于包围球与包围盒相结合的更准确的可见性剔除方法和连续无缝化的细节过渡方案，有效地避免了绘制失真和模型的闪烁现象。

4.基于可编程GPU技术实现了地形渲染的加速和增强，即通过使用着色器语言控制图形绘制流水线的各个阶段，将高程数据的顶点信息和索引信息分别生成两张纹理存入显存，供整个地形绘制调用；利用曲面细分和分形技术在几何阶段对顶点插值和偏移，生成过程化细节，弥补分辨率不足时地形网格棱角化的现象。

附图说明

图1是UTM投影图；

图2是多分辨率金字塔基本结构图；

图3是全球多分辨率地形分层分块方案图；

图4是滤除噪声的流程图；

图5是滤除噪声的效果图；

图6是使用不同的步长值从相同的顶点缓冲获取的不同细节层次图；

图7是动态节点的生成过程图；

图8是细分过程的示意图；

图9是像素值的确定图；

图10是噪声图的合成图；

图11是两种算法生成的噪声图；

图12是细节合成效果图；

图13为两组数据的高程图元数据；

图14为两个场景的鸟瞰图；

图15为本发明实现的LOD无缝连续绘制效果图；

图16为渲染的实时细节增强效果图；

图17为利用GPU监测工具NVIDIA nTune和GPU-Z对系统的运行性能进行监测获得的记录图；

图18为两组实验按照固定路线漫游时的三角形数目和相应的帧速变化图；

图19为几何误差转换至屏幕误差图；

图20为透视投影中的视景体图；

图21为视椎裁剪后的结果图；

图22为6级LOD的细节分配方案图；

图23为网格的几何过渡图；

图24为多级几何过渡图。

具体实施方式

1.研究了对地形数据进行组织管理的预处理流程，从诸多的表示模型和高程所采用的投影坐标系开始，依次讨论了投影变换、噪点消除和插值采样的实现方法，为构建多分辨率地形数据金字塔做准备。

2.研究了地球空间剖分的方案和多分辨率金字塔分层分块的数据组织方式，对全球多分辨率地形数据按照本文设计的文件格式进行了组织存储；使用四叉树结构对各个数据节点构建索引，并提出了四叉树的更新插入快速算法。

3.研究了基于GPU的三维地形绘制增强技术。分析了图形绘制流水线的各个阶段的操作与分工，讨论了地形数据与GPU编程技术相结合的绘制算法，充分利用了GPU的并行计算能力，实现缓存复用和实时细节增强，大大降低CPU和内存的负载，提高了绘制性能和显示效果。

4.基于本文的研究内容，给出了全球地形数据管理与显示系统的概要设计模型，划分了各个模块的功能，并讨论系统的工作流程，用实验加以验证。

本发明绘制出三维地形模型的多分辨率金字塔模型构建：

①地形数据的预处理流程：通过Ⅰ、投影坐标的变换，Ⅱ、影像噪点的消除和Ⅲ、采样插值构建多分辨率金字塔统一了数据格式；

Ⅰ、投影坐标的变换：对UTM投影的反解：

实时的三维显示系统会采用一种笛卡尔右手坐标系来统一在世界空间内的模型。但考虑到建立全球地理参考系统，这种坐标系比较接近于以地球质心为原点的空间直角坐标系，所有的位置坐标均是以米为单位，用(x，y，z)表示相对于质心的偏移。但这种坐标的值通常很大，大多超过了单浮点型数据的表示范围，单浮点型数据的格式中仅能保存23个比特位的整数部分(最高的一位作为正负标识符)，转化成十进制也就是2²³-1＝8388607，而小数部分最多只有7位。在WGS84参照系中，地球的赤道半径就达到了6378137m，这样的话使用单浮点型格式，最高空间分辨率也不会超过0.8m(6378137/8388607≈0.8)。在这样的阈值下，各种各样的浮点数据近似解所造成的失真就会发生，如顶点重叠和视点抖动等。

另外一个因素是，空间直角坐标系无法直观地表示出制图的精度信息，相反，很多地理或是投影坐标系却可以表示更多的地理信息，地理坐标系统与其所采用的参考椭球相关，如经纬度地理坐标；投影坐标系统则将地球模型投影到一些简单的展开曲面上，Lambert Conformal Conic(LCC)采用的是正轴等角圆锥投影，Universal TransverseMercator(UTM)采用的是横轴等角割圆柱投影。因而将地球空间投影在直角坐标系中的方法并不可取，这同样需要投影坐标系与地理坐标系之间的正反变换来确定点的位置。

作为地形数据源的DEM采用的是WGS84作为参考椭球体，所采用的坐标格式有两种，一种是UTM坐标，单位是长整型的米或千米，一种是WGS84经纬度坐标，单位是双精度类型的角度。为了统一数据格式，根据人们的实际需求，组织DEM数据时，采用经纬度坐标作为我们的度量单位。因此，需要将以UTM坐标为单位的DEM数据转换为经纬度坐标。

对UTM坐标进行反解，即是由其投影方程：

其中，x代表纬度的投影值，y代表经度的投影值；

设中央经线所在的子午圈经度值为L₀，对应每一条纬圈，除了赤道投影为直线外，其余均向内凸起，对X＝x直线与中央经线的交点称之为底点F，其经纬度坐标为(B_F，L₀)，这条直线上点的经纬度坐标以F点为原点的表达式设为：

进而B＝B_F+B′，L＝L₀+L′，原式就转化为了求解B_F，B′和L′三个参数的过程；由于直线上的点相对于底点为对称投影，且同一条直线上x的值不变，则B′必为关于y的偶函数，L′必为关于y的奇函数，按照傅里叶函数展开可得：

其中，n_i(i＝0，1，2，...)是关于x的函数，而也就是

由可得到：

n₀的求解可由高斯投影的算法得出，进而可依次求出公式中的各项系数，最终可得：

其中各个参数的表示如下：

Ⅱ、影像噪点的消除：

通过分析高程图中的噪声特征，可以发现，高程的噪声主要是白噪声，而且噪声面积较大，在影像中呈块状，而不是单纯的盐噪声，噪声模型符合加性噪声的特征：

其中，I(x，y)为位于影像(x，y)处的灰度值，O(x，y)为未受噪声污染的灰度值，n(x，y)是独立于O(x，y)的随机噪声值，p为受到噪声污染的概率；

因此在滤波时，不再考虑极小值作为噪声点，而且将滤波窗口的初始尺寸设定的较大，以避免漏检块状噪声。

在处理之前，需要首先明确一下几个定义：

拉格朗日的插值方法为：

其中，y_i表示x_i出的函数值；此处，插值算法是针对一维数据进行处理，当拓展到二维影像上时，我们取插值像素点(x，y)水平、垂直、两个对角线方向共四个方向上的临近像素点进行插值，并取其平均值作为最后的插值结果；

算法的实验效果如图5所示，可以看到噪点被全部检测出并成功滤除，尽管噪点周围的高程受其影响被提升了一定的高度，但相对于中值滤波后整个地形被平滑的效果该方法还是可以接受的。

对影像进行滤波的步骤流程如下：

(ⅰ)、设置初始的滤波窗口大小basicWndSize为11×11，最大的窗口大小maxWndSize为19×19，设置噪点计数器noiseCount＝0，以及表征像素点是否是噪点的标识符isNoise＝false，计算出整幅图像的均匀度image_HSV，然后对整幅图像的像素进行遍历；

(ⅲ)、判断当前窗口中心像素点(x，y)的灰度值I(x，y)是否等于窗口内像素极大值iMax，或者窗口均匀度wnd_HSV是否大于影像均匀度image_HSV；若为是，则继续判断当前窗口内极大值像素的个数nMax是否小于sumPixel/3，若为是，则标记中心像素点为噪声点isNoise＝TRUE，noiseCount++；若为否，则认为当前的极值密度过大，进而影响到对噪声的判断，此时将滤波窗口扩大basicWndSize+＝2，跳转到步骤2；如果滤波窗口已经达到了最大值，即basicWndSize＝＝maxWndSize，则标记中心像素点为噪声点isNoise＝TRUE，noiseCount++；否则，对该店进行中值滤波；若为否，则判定当前点不是影像噪声点，则跳出循环，将滑动窗口移至下一个像素，直至对整幅图像完成遍历；

(ⅳ)、对判定为噪声点的四个方向进行拉格朗日插值；

Ⅲ、采样插值：即DEM数据的插值采用双线性插值方法；取其平均值作为噪点的像素值，将该点的标识符isNoise设为false，同时noiseCount--。

②构建多分辨率金字塔：

三维地形在进行实时绘制时，为了保证视觉效果的连续性，需要对原始地形网格采样或插值，形成不同分辨率表达的网格，根据算法预先设定的误差度量控制相应分辨率网格的调用或移除，但这种多分辨率网格难以通过实时计算获得并应用于图形硬件的渲染上，网格尺寸较小时尚不足以引起帧速的下降，但当视点收容范围较大时对帧速的影响就显得不容忽视。

为了解决此问题，产生了多分辨率金字塔的概念，金字塔结构是一种多分辨率层次模型，从模型的原理上讲，金字塔是一种连续分辨率模型，但在构建时难以实现分辨率的连续变化，算法的复杂度和空间的冗余会大大提高，因此，金字塔构建时，多采用倍率方法构建，即对同一地区，采用一系列的网格表示，上下相邻网格的采样间隔相差相同的比例系数m(通常取m＝2)，这样从金字塔自底向上，分辨率逐层降低，网格尺寸也相应减小。设原始网格的尺寸为S₀，分辨率为R₀，所构建的金字塔共有l层(金字塔顶端为第0层)，那么第k层的网格尺寸s和分辨率r可以表示为：

对平面投影后的地球按照等经纬度进行分块，同时按照自顶向下的方式，构建不同层次级别的金字塔层，见图3，构建的方案如下：

(Ⅳ)每个分块中包含了对应地理位置的DEM数据，其分辨率尺寸必须是奇数，即(2n+1)×(2n+1)；以保证与其邻块有一个像素的重叠率，这样不仅便于避免T型裂缝的出现，对高程数据插值时，可以使相邻的地形块连接处高度一致，消除了棱台现象的发生；

(Ⅴ)金字塔的底层为最高分辨率，直接从原始数据中获得，其所在的LOD层级由其自身分辨率与所属范围决定；每个地形块对应的纹理数据和法向数据分辨率必须要高于DEM的分辨率，以保证插值平滑，避免渲染失真；假设某一层DEM分辨率为r，这里把它们的分辨率定义为r′＝n×(r-1)，n＞1，n∈N；得到了全球规模的金字塔模型。

其属性参数如下表所示：

表金字塔模型的分层分块结果

其中，R表示WGS84所定义的地球长半径：6378137.0m。对全球数据分块的依据是经纬度的跨度差，而表中第4、5列的单位却是以距离表示，以更好地符合人们的使用习惯，它们的值取自于所在块跨越的最长纬度距离及其与像素分辨率的比值。

本发明基于GPU的流水线绘制增强技术：动态图像处理：

a、节点的数据结构与组织：

将世界分割为许多的小区域，每个都以长方体表示，称之为区块(sector)，将这些区块形成一个层次树，每个区块都由四个更小的子区块构成。进行地景渲染时，首先判断哪些区块是可见的，如果父区块不可见，那么孩子也不可见；反之则孩子节点都可见。而当父区块与视椎体相交时，就选择它的子区块，将它继续当作父区块来测试其可见性。直到所有可见块被找到或完成了所有节点的遍历。

取一张2D的高程图为例，首先对其使用Geomipmap进行滤波，构建一张包含了L层的mipmap金字塔。在这里，我们将每层的数据都以适当的尺寸保存为2D纹理，而不是认为的通过顶点数组将之线性化为1D的缓冲。金字塔的每一层各包含了n×n的采样点，我们为之存储了一张单通道的高程图，一张三通道的顶点法向图和三通道的纹理图。

这样，我们可以在着色器中定义如下的节点数据结构：

让法线图与纹理的分辨率为高程的两倍，因为每个顶点只有一条法线的效果很模糊，同时在加入光照后也会很粗糙，而两倍分辨率的法线图可以通过插值创造更加柔和的模型。低分辨率的纹理限制视点的范围，使各个表面类型产生大面积的模糊，故纹理分辨率也不应小于高度图的分辨率。

建立两个类，一个是cTerrain，它代表整个地形，来帮助裁剪出实际需要的部分地形，而这个地形被划入网格对应的区块。这些区块使用第二个类表示，cTerrainSection，此类中存储了每个区块的高程及法线。例如，一个256×256个顶点的地形，我们将创建cTerrain来存储整个数据集。cTerrain对象进一步将地形细分为32×32的顶点区块，创建总共64个cTerrainSection对象来表示每个区块。既然所有的cTerrainSection对象都是等尺寸，等数量顶点的，我们就可以建立一个单独的索引缓存对象并使用它来渲染地形的各个区块，针对不同的层级，只需将相应的步长值向右移一位，如图6所示。

构建索引缓冲的代码如下所示：

bool cTerrain::buildIndexBuffer(

int gridVerts,//网格边上的顶点个数

int step)//索引缓冲使用的步长值

{int totalStrips＝gridVerts-1；//条带的个数

int indice_tripwise＝gridVerts<<1；//每个条带占两行顶点

//由于0三角的存在，顶点个数增加

int total_indice＝(indice_stripwise*totalStirps)+(totalStrips<<1)-2；

int*index＝new int[total_indice]；

int vert,start_vert＝0；

for(int j＝0；j<total_strips；++j)

{//生成一行的索引缓冲

vert＝start_vert；

for(int k＝0；k<xVerts；++k)

{*(index++)＝vert；

*(index++)＝vert+gridVerts*step；

vert+＝xStep；

}

start_vert+＝lineStep；

if(j+1<total_strips)

{//在每一行的末尾增加0三角

*(index++)＝(vert-xStep)+lineStep；

*(index++)＝start_vert；

}}}

构建顶点缓冲的具体的做法，即将每个块的位置和纹理坐标(x，y)收缩在(0，1)的范围内。在渲染时，我们将坐标放大至需要的比例，加以合适的偏移值即可移动到对应的世界位置。而各个地形块使用一系列相似的xy位置和uv纹理数据，我们将它们一次性存储与cTerrain对象内渲染每个cTerrainSection对象时再引用它们。这样我们构建了两组顶点数据流：存储于cTerrain的共享顶点数据和各个cTerrainSection的独立顶点数据，大大减少了其占用的内存空间。

构建顶点缓冲代码如下：

bool cTerrain::create(

cTexture*heightMap,

const cRect3d&worldExtents,

int shift)

{

//设置放缩比例

m_tableWidth＝(uint16)heightMap->width()；

m_tableHeight＝(uint16)heightMap->height()；

m_mapScale.x＝m_worldSize.x/m_tableWidth；

m_mapScale.y＝m_worldSize.y/m_tableHeight；

m_mapScale.z＝m_worldSize.z/255.0f；

//设置初始偏移

m_sectorCountX＝m_tableWidth>>m_sectorShift；

m_sectorCountY＝m_tableHeight>>m_sectorShift；

m_sectorSize.set(m_worldSize.x/m_sectorCountX,

m_worldSize.y/m_sectorCountY)；

cVector2vert(0.0f,0.0f)；

sLocalVertex*pVerts＝

new sLocalVertex[m_sectorVerts*m_sectorVerts]；

//将共享的顶点与纹理坐标写入流中

for(int y＝0；y<m_sectorVerts；++y)

{

vert.set(0.0f,y*cellSize.y)；

for(int x＝0；x<m_sectorVerts；++x)

{

pVerts[(y*m_sectorVerts)+x].xyPosition＝vert；

pVerts[(y*m_sectorVerts)+x].localUV.set(

(float)x/(float)(m_sectorVerts-1),

(float)y/(float)(m_sectorVerts-1))；

vert.x+＝cellSize.x；

}}}

c、三维地形模型的实时绘制：分为两类，一种是静态节点，另外一种是动态节点；

Ⅰ、大规模地形数据渲染：即在预处理过程中通过采样获取同一区域不同分辨率节点，在渲染时，通过将对应分辨率的节点调入显存中直接进行渲染；

采用Chunked LOD的四叉树结构组织地形数据，同时结合Geo-Clipmap中的嵌套网格模式进行绘制，并使用一种连续的几何过渡来调整LOD变化时网格的分裂过程，避免突变和裂缝的产生。

①环状取址的位运算：

视点移动后的像素范围变成了(x_min+x_diff，y_min+y_diff)和(x_max+x_diff，y_max+y_diff)，当更新节点队列时，从节点到队列映射需要满足以下三个条件：

(1)同一个元素i始终对应队列中的同一个元素f(i)，即从i到f(i)满足1:1的映射关系；

(2)对任意元素x，一定存在f(x)与之相对应，保证映射的连续性；

(3)若那么保证映射的无冗余性。

第一个条件保证了增量更新的可能性，一旦节点加入到了队列中时，如果节点在视点移动后仍保留在视场内，那么就不需要将之移出队列；第二和第三个条件防止数据被改写，使队列保持最小空间。实现环状取址，利用模运算实现上述映射。

其实现代码如下所示：

其中，(i+x_min+x_diff，i+y_min+y_diff)表示节点的实际坐标，(x₀，y₀)表示节点在队列中的坐标。这样一来，无论是视点向任何方向运动时，均可以通过一个函数读取到对应位置的坐标，省去了在循环中多次判断位移正负方向的繁琐步骤。

当视点位置向右上方移动一个节点的宽度时，节点ABCDEFG进入active region，与此对应的左下角的“L”形结点会被移除，在节点缓存中，底边节点和侧边节点分别被取代。当继续向右上移动一个节点的宽度时，节点HIJKLMN继续进入active region，除了第二行与第二列原有的节点被移除以外，还应注意到节点A和G也同时被移除了，此时的队列数组按照底纹不同，被划分成了四个区域，除了右上角的原有节点保持不变外，EFMN为一个区域，JKDL为一个区域，HIBC也组成了一个区域，环状取址的方式使“L”形更新区域在队列中变成了“+”形，因而在读取更新节点时还可以按照这三个区域划分进行读取。

②误差度量的确定：

采用Chunked LOD的误差度量计算方式，即从叶子节点开始，沿着四叉树向上逐层计算节点的误差，保证父节点的误差度量始终大于子节点，对于每个节点的误差度量∈_m，取值于max{ε₁，ε₂，…，ε_n}，其中n为节点所包含的网格内被移除顶点的总数，ε_i表示简化过程中移除一个顶点后地形产生的高度差值。而为了保证误差的嵌套，一个节点的误差度量ε_c为其四个子节点的ε_a的最大值加上其本身的误差度量ε_■；ε_c的计算法如下：

误差ε是基于三维物空间测量的几何误差，当利用透视投影将几何误差投影在二维的屏幕空间时，实时的LOD选择可以更好地利用二维的屏幕空间误差阈值。当确定视点参数后，根据当前的几何误差ε确定对应的屏幕误差ρ，如图19所示，ω表示视见体的宽度，d表示视点到当前层的距离，θ为视场角大小，x为屏幕分辨(以像素为单位)。由相似三角形可得出：

③可见性剔除：

地形网格的选取是一个重复迭代的过程，会不断地重复对每个节点网格计算距离与误差，这样就涉及到大量的求交运算。为了减少待处理的网格数量，采用包围盒与包围球相结合的判断方式提高运算速度。一个物体的包围体是能把这个物体包裹在内的一种形状，通过快速地排除与视椎体不相交的包围体，即可减少很多不必要的运算。如果视椎不与包围体相交，那么它必然不会与包围体内部的物体相交，而选取形状简单的包围体要比判断其内部外形复杂的物体要节省很多时间。

视点较近时的剔除方法：

当视点距离地形较近时，采用包围盒与视椎体进行相交测试，即通过包围盒的特征点与视椎体的6个平面依次求交，以确定这些特征点是否属于平面正向的那半个空间或者与平面相交；这就意味着视椎体空间是由这些平面的正半空间组成的。因此，如果一个矩形存在于任何视锥平面的负半空间的话，那么这个物体就在视锥之外了。视椎的6个平面由方程f(x，y，z)＝ax+by+cz+dw＝0决定，其中向量(a，b，c)表示该平面的法向量，数值dw表示平面与原点的距离；包围盒的9个特征点由8个角点和一个中心点组成，通过比较包围盒的特征点即可判断地形网格的可见性，取size等于包围盒的斜对角线长的一半；

对于空间内任意一点P(x₀，y₀，z₀)与平面f(x，y，z)＝ax+by+cz+dw＝0的关系，计算公式：

结果的绝对值即表示点到平面的垂直距离，比较绝对值与size的大小，就能得到节点是否在视椎内部。

视点较远时的剔除方法：

f(x，y，z)＝m(x-x_i)+n(y-y_i)+p(z-Z_i)＝0 (3.9)

mx+ny+pz＝0 (3.10)

计算出所有8个顶点的p_i′后，取8个坐标点中的极大极小值min X，max X，min Y，max Y，min Z，max Z；令

距离系数a的值为：

④缝隙消除与几何过渡：

提出了一种新的LOD层级转换的控制技术，同时删去在不同分辨率网格之间的退化三角形插入操作，避免不必要的渲染负担和由此引发的渲染失真，以实现无缝化的LOD过渡。

为了控制在一定的距离下选取适当的细节层次，在进行节点筛选前，首先计算每个LOD层级所覆盖的地理范围，以便渲染时在视窗平面上产生颗粒度近似一致的三角形数量；由于四叉树结构的定义，每层细节与它上一级的顶点数目之比固定为4，因而在一维方向上的距离倍数也就为2，因此不同分辨率的过渡系数设定在2左右以保证三角形数目近似一致。图22即为不同的LOD层级所覆盖的距离范围，随着分辨率降低，每层地形块所占的地理范围是上一级的2倍，而最后一层(距视点最远，细节层次最低)即代表整个视点的可视范围。

为了保证相邻LOD层级之间能够平滑过渡，同时尽可能地保留地形精度，我们采用与Geo-Clipmap算法一致的过渡方案，即在低分辨率网格与高分辨率网格相邻的区域中划分出一部分作为过渡区域，由低向高过渡，实验中选取每个LOD范围的33％作为过渡区域。这样做的副作用是会产生更多的三角形，但不同于退化三角形的生成方法，由顶点位移方法生成的三角形在复杂度上要简单的多。

所有过渡区域的顶点均是按照自身唯一的误差度量进行变换，即逐点进行过渡，每个顶点均支持两种变换方式：由本身向更高分辨率的过渡和由本身向更低分辨率的过渡。取第二种变换方式为例，过渡操作即是每个3×3的顶点网格，由8个三角形组成，通过逐步地移动4个边中点和1个中心点到左下角点的位置，使之简化为2×2的顶点网格，由2个三角形组成。整个过程看起来，如同右上角的两个三角形逐渐扩大，而剩余的6个三角形一点点向左下方向压缩成为退化三角形并最终不再栅格化为像素的过程。这种方法能够有效地消除相邻LOD层级之间的缝隙，并产生平滑的过渡，不会出现T型连接。

首先，计算顶点到视点的近似距离来决定位移系数的大小，用来计算距离的顶点位置根据网格在世界坐标系中的位置比例近似获取，因为网格在坐标内连续排布，则获取的顶点位置也应该固定，相邻网格的边缘顶点必须一致，这样才可以保证网格之间不会产生缝隙。在实验中，x和y分别取自经纬度值，而高程值则是由高程图连续采样得出的。由图22，设每层LOD的起始距离为start_i，i表示当前LOD所在的层级，其所覆盖的距离范围为range，而视点到当前LOD任意一个顶点的距离为dist，则当前的过渡区域[mStart,mEnd]以及该顶点的过渡系数k为：

其中，clamp(x，a，b)为钳位函数，即

这样一来，k的大小即被限定在了0(完全没有变换，停留在低分辨率位置)和1(完全变换，三角形数变为原来的两倍)之间，而随着k的值从0到1逐渐变大，网格也实现了由低分辨率向高分辨率的过渡。在获得过渡系数后，定义顶点所在网格尺寸为n×n，所在的四叉树深度为l，顶点在网格内归一化的坐标为(x，y)，则网格尺寸相对于世界坐标下的实际距离放缩值scale以及过渡后的坐标(x′，y′)由如下公式求得：

其中，frac(x)表示取x的小数部分；

最终，顶点坐标会从(x，y)变换到(x/2，y/2)的位置，它对应的高程通过双线性插值获得；当所有顶点过渡完成后，当前的网格就会与它下一层的高分辨率网格完全重合，而期间的过程不会产生任何缝隙。图23为不同的过渡系数对应的网格形状。

事实上，由于视点是连续变化的，而整个视场中，所选取的节点分辨率随着视景的加深逐渐降低，因而每一帧都会发生不止一层的几何过渡，如图24所示，多层过渡带能够紧密地连接相邻分辨率网格，避免裂缝的产生。实验表明，在进行地形漫游时，网格间平滑连续的过渡方法有效地消除了节点细分时造成的“突跳感”，帧速也稳定在60帧以上，不会产生视觉上的停滞。

Ⅱ、动态节点，它是在四叉树结构延伸至叶子节点时，仍不能满足视点阈值的需求，由GPU实时计算生成并加入四叉树中的节点。

动态节点的生成：

地形模型的细分：细分是在四边形网格的原始顶点上根据细分掩膜插入面点和边点，并与原始顶点连接成新的四边形，生成的非边界顶点度数为4，边界顶点度数为3，都是正则点；如图8所示，正则点细分过程：对四边形V₅V₆V₇V₈进行细分时，会在其四条边上形成四个边点E₀E₁E₂E₃，同时还有一个面点F，使原来的面片分裂成四个子面片；E₀，E₁，E₂，E₃的值如下公式：

面点F的值由公式4.2得出：

其中，

这种细分算法只能应用于正则点所在规则格网上，对于奇异点不予考虑。由于Chunked LOD算法的特点，细化过程是在每个块内独立进行的，对于邻块内的网格无从知晓，当掩膜移动至块的边缘时，传统的边外推和对角点替换都会产生假边缘，为了保证视觉上的真实感，我们在块结构中再存储上邻块的边采样点，只有在对边进行细化时才使用其作为控制点，块内的细分则对其不予考虑，以增加少量数据为代价换取边缘的连续性。

随机分形：用变量(x，y)表示曲面上点的坐标，X(x，y)表示该点在曲面上的高；则FBM曲面定义为在某概率空间上，指数为H(0＜H＜1)的一个随机过程

随机变量X：R→R²满足：

使用分形技术绘制的三维地形能够产生无限具有自相似性的逼真场景，易于进行平滑和插值。这里采用两种分形算法来生成不同样式的噪声图，分别是Diamond-Square算法和Perlin噪声。

1.Diamond-Square算法

算法的实质是二维平面上的随机中点位移算法。它由种子点组成的正方形开始，通过若干次中点位移法迭代后，不断细分种子正方形以获得逼真的三维地形。在迭代过程中每一次细分都会执行两步操作，Diamond步和Square步。

Diamond步用于生成正方形的中心，取四个顶点计算平均值，再加上一个随机量，作为中心点的值。Diamond步生成了四个三棱锥。Square步则是取棱锥的四个顶点，用其平均值加上随机变量(与Diamond步的随机变量服从相同的概率分布)作为原正方形的四条边中点，每一次细分的结果同图8一致。

在过程的各个阶段，随机值总是生成在固定范围内，我们称当前范围为delta，总范围是[-delta,delta]。在下一次的细分时，delta会乘上一个放缩值，称之为粗糙度(roughness)，介于0到1之间，用来减少delta的取值范围。高粗糙值会制造出更加无序的地形，而低粗糙度会使地形越来越平坦。

Diamond-Square算法实现简单，运算速度较快，被广泛地采用。但由于相邻正方形间没有信息传递，因而导致噪声图中产生“折痕问题”，即图像中存在不自然的线性轨迹，并且不能通过局部平滑移除；另外，其迭代生成方式计算量较大，耗费过多时间而不够灵活。

2.Perlin噪声

Perlin噪声是一个生成随机噪声的函数，Perlin噪声通过对多重分形的合成获取最终结果，即利用有限带宽频率的逐次增加合成，而各个频段的计算在合成过程中相互独立，因此它对地形分形的频率、维数以及其它相应统计特性具有较强的控制能力。

对于一幅噪声图像，首先给整数坐标点赋予随机方向的法向量与灰度值，除此以外任意一点的像素值由该点所在最邻近的单元格内四个顶点的法向量，和连接该点与四个顶点的方向向量来决定。

设为图中所求点的向量，n为噪声的叠加次数，一般取n∈[3，12]。设噪声函数波长为λ，则频率ω＝λ^-0.5i，该点的灰度值为：

其中，offset为偏移量取-0.1，是一个噪声函数。下图显示了两个独立噪声加和后的样子。

采取预处理的方法生成多种类型的噪声纹理，避免实时计算增加GPU的额外开支，在光栅化阶段中，像素着色器读取一张包含噪声的纹理图与高程图进行融合，这样的执行效果更加快速。同时噪声图进行不同比例地放缩与叠加，就可以扩展出无限大的噪声纹理，以此消除绘制出的周期性和模式化效果。图11显示的是两种噪声算法所生成的噪声纹理图，图12则是将两种噪声分别与高程叠加后所产生的效果。由图可以看出，二者可以实时地增强过程化细节，模拟出不同的地形效果，Diamond-Square噪声更加接近沙地效果，而Perlin噪声则接近岩石效果。

实现效果：

两组实验数据，分别构建了8级LOD的世界高程与12级LOD的Puget Sound的三维地形，其中8级世界高程的分辨率为152m，覆盖地理范围整个南北纬83°以内的整个地球空间，构建而成的原始网格由268,435,456个三角形组成，而12级的Puget Sound的高程分辨率为10m，覆盖的地理范围是163km×163km，构建而成的原始网格由463,601,664个三角形组成。由于高程信息的数据格式并没有采用压缩，12级的LOD文件大小达到了2.58GB，其对应的法向图则接近有4GB。

下表为开启细节增强前后系统性能的变化。数据显示，随着视点的降低，LOD层级也逐渐深入，三角形数目不断增长，但同时与视见体相交的块也减少了，所以第8级LOD的三角形数目在高度较低的情况下突然减少。第9级和第10级是动态节点加入后的结果，三角形数目开始增加，但载入速度也有所减慢。同样，前8级是由预处理过程生成的静态节点，浏览时仅需从存储中读取对应的地形块，帧速较为稳定，而动态节点的插入给实时渲染带来一定负担，帧速的变化取决于地形复杂度和三角形数目。

表静态节点与动态节点的运行时间比较

视点高度(m)	视点所在LOD层级	三角形个数	载入时间(s)	漫游帧速(s^-1)
					1846.2	5	45056	8.96	145.22
1282.8	6	40960	7.82	147.63
					1015.5	7	32768	7.46	146.97
882.2	8	6144	5.88	159.46
					424.7	9	30720	7.52	133.23
301.1	10	26624	7.33	140.60

由图17两组监测数据可以看出，GPU承担了大部分的运算负载，CPU和内存的占用率一直处在较低的水平，由于两组数据所使用的四叉树深度不同，视场范围内所容纳的节点个数也不尽相同，这是造成显存占用不一致的主要原因，但相比原始文件GB级的大小，本系统已经实现了利用少量缓存浏览大规模地形的功能。

统计所得的数据，系统在运行时平均每一帧所绘制的三角形数目达到265.6k，而帧速则接近194s^-1，由此可知本文利用GPU对绘制的效率提升效果明显，而由图18看出，帧速与三角形数目并不是严格的负相关关系，随着四叉树的深度加深，地形的细节越来越丰富，会导致帧速的明显下降，即帧速受地形的复杂程度影响也较大，8级的世界高程多为平坦的地形，因此其帧速受三角形数目的影像较小，12级的Puget Sound则波动较大。下表列出了一些算法的帧速比较，但由于不同的算法采用的硬件平台不同，而本文的算法对GPU的依赖性较强，因此难以进行横向比较，在下一步的研究中应尽量降低算法所造成的绘制硬件负载。

表不同算法的绘制效率比较

算法	平均帧速(s^-1)
		经典算法	10-20
ChunkedLOD	21-50
		Geo-Clipmap(2004)	95-120
Geo-Clipmap(2005)	130-185
		CDLOD	56-1800
本文算法	82-280

本发明进行地形多分辨率显示验证系统的设计与实现。首先分别从框架结构和各个模块的功能两个方面介绍了系统的设计方案，然后用两组实验数据验证了系统的主要功能，并用图表加以展示，最后经过专业的检测软件进行了性能测试，证明该系统满足了实时和GPU友好的绘制要求。

Claims

1.一种结合GPU技术实时流畅绘制出三维地形模型的方法，其特征在于：

三维地形模型的多分辨率金字塔模型构建：

Ⅰ、影像噪点的消除：

噪声模型符合加性噪声的特征：

对影像进行滤波的步骤流程如下：

<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>d</mi> <mrow> <mo>(</mo> <msub> <mi>I</mi> <mrow> <mi>M</mi> <mo>&times;</mo> <mi>N</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>M</mi> <mo>&times;</mo> <mi>N</mi> </mrow> </mfrac> <munder> <mo>&Sigma;</mo> <mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>&Element;</mo> <msub> <mi>I</mi> <mrow> <mi>M</mi> <mo>&times;</mo> <mi>N</mi> </mrow> </msub> </mrow> </munder> <mfrac> <mrow> <mo>|</mo> <mi>f</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>a</mi> <mi>v</mi> <mi>e</mi> <mi>r</mi> <mi>a</mi> <mi>g</mi> <mi>e</mi> <mrow> <mo>(</mo> <msub> <mi>I</mi> <mrow> <mi>M</mi> <mo>&times;</mo> <mi>N</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mrow> <mi>a</mi> <mi>v</mi> <mi>e</mi> <mi>r</mi> <mi>a</mi> <mi>g</mi> <mi>e</mi> <mrow> <mo>(</mo> <msub> <mi>I</mi> <mrow> <mi>M</mi> <mo>&times;</mo> <mi>N</mi> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>d</mi> <mrow> <mo>(</mo> <msub> <mi>W</mi> <mrow> <mi>r</mi> <mo>&times;</mo> <mi>r</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>r</mi> <mo>&times;</mo> <mi>r</mi> </mrow> </mfrac> <munder> <mo>&Sigma;</mo> <mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>&Element;</mo> <msub> <mi>W</mi> <mrow> <mi>r</mi> <mo>&times;</mo> <mi>r</mi> </mrow> </msub> </mrow> </munder> <mfrac> <mrow> <mo>|</mo> <mi>f</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>a</mi> <mi>v</mi> <mi>e</mi> <mi>r</mi> <mi>a</mi> <mi>g</mi> <mi>e</mi> <mrow> <mo>(</mo> <msub> <mi>W</mi> <mrow> <mi>r</mi> <mo>&times;</mo> <mi>r</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mrow> <mi>a</mi> <mi>v</mi> <mi>e</mi> <mi>r</mi> <mi>a</mi> <mi>g</mi> <mi>e</mi> <mrow> <mo>(</mo> <msub> <mi>W</mi> <mrow> <mi>r</mi> <mo>&times;</mo> <mi>r</mi> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2.15</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <mi>a</mi> <mi>v</mi> <mi>e</mi> <mi>r</mi> <mi>a</mi> <mi>g</mi> <mi>e</mi> <mrow> <mo>(</mo> <msub> <mi>I</mi> <mrow> <mi>M</mi> <mo>&times;</mo> <mi>N</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>M</mi> <mo>&times;</mo> <mi>N</mi> </mrow> </mfrac> <munderover> <mo>&Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>M</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>N</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <mi>f</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2.13</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <mi>a</mi> <mi>v</mi> <mi>e</mi> <mi>r</mi> <mi>a</mi> <mi>g</mi> <mi>e</mi> <mrow> <mo>(</mo> <msub> <mi>W</mi> <mrow> <mi>r</mi> <mo>&times;</mo> <mi>r</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>r</mi> <mo>&times;</mo> <mi>r</mi> </mrow> </mfrac> <munderover> <mi>&Sigma;</mi> <mrow> <mi>m</mi> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <mi>r</mi> <mo>-</mo> <mn>1</mn> </mrow> <mn>2</mn> </mfrac> </mrow> <mfrac> <mrow> <mi>r</mi> <mo>-</mo> <mn>1</mn> </mrow> <mn>2</mn> </mfrac> </munderover> <munderover> <mi>&Sigma;</mi> <mrow> <mi>n</mi> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <mi>r</mi> <mo>-</mo> <mn>1</mn> </mrow> <mn>2</mn> </mfrac> </mrow> <mfrac> <mrow> <mi>r</mi> <mo>-</mo> <mn>1</mn> </mrow> <mn>2</mn> </mfrac> </munderover> <mi>f</mi> <mrow> <mo>(</mo> <mrow> <mi>i</mi> <mo>+</mo> <mi>m</mi> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2.14</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>L</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <mi>x</mi> </munderover> <msub> <mi>y</mi> <mi>i</mi> </msub> <msub> <mi>l</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>l</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <munder> <mi>&Pi;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>0</mn> </mrow> </munder> <mrow> <mi>j</mi> <mo>&NotEqual;</mo> <mi>i</mi> </mrow> <mi>n</mi> </munderover> <mfrac> <mrow> <mi>x</mi> <mo>-</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> </mrow> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> </mrow> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2.16</mn> <mo>)</mo> </mrow> </mrow>

其中，y_i表示x_i处的函数值；此处，插值算法是针对一维数据进行处理，当拓展到二维影像上时，取插值像素点(x，y)水平、垂直、两个对角线方向共四个方向上的临近像素点进行插值，并取其平均值作为最后的插值结果；

②通过上述消除噪点后的数据构建多分辨率金字塔：

<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>s</mi> <mo>=</mo> <msub> <mi>S</mi> <mn>0</mn> </msub> <mo>&times;</mo> <msup> <mi>m</mi> <mrow> <mi>k</mi> <mo>-</mo> <mi>l</mi> </mrow> </msup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>r</mi> <mo>=</mo> <msub> <mi>R</mi> <mn>0</mn> </msub> <mo>&times;</mo> <msup> <mi>m</mi> <mrow> <mi>k</mi> <mo>-</mo> <mi>l</mi> </mrow> </msup> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2.7</mn> <mo>)</mo> </mrow> </mrow>

其中，m＞1；

(Ⅳ)每个分块中包含了对应地理位置的DEM数据，其分辨率尺寸必须是奇数，即(2n+1)×(2n+1)；

(Ⅴ)金字塔的底层为最高分辨率，直接从原始数据中获得，其所在的LOD层级由其自身分辨率与所属范围决定；每个地形块对应的纹理数据和法向数据分辨率必须要高于DEM的分辨率，以保证插值平滑，避免渲染失真；假设某一层DEM分辨率为r，这里把它们的分辨率定义为r′＝n×(r-1)，n＞l，n∈N；得到了全球规模的金字塔模型；

基于GPU的流水线绘制增强技术，对上述全球规模的金字塔模型的数据进行绘制：

a、节点的数据结构与组织：

采用ChunkedLOD算法构建了块状四叉树结构，即不再以单个顶点或三角形为处理单元，而是以包括纹理贴图在内的网格模型，将三角形打包成条带一次性完成渲染，充分利用GPU的批处理能力；

①环状取址的位运算：

移出视场的节点被第一次进入视场的节点代替，其余节点不动，这样就完成了增量更新，对于视场范围内任意一点(x，y)，其在节点数组中存储的位置为(xmodn，ymodn)，其中“mod”表示求模运算；

②误差度量的确定：

采用ChunkedLOD的误差度量计算方式，即从叶子节点开始，沿着四叉树向上逐层计算节点的误差，一个节点的误差度量ε_c为其四个子节点的ε_ci的最大值加上其本身的误差度量ε_m；ε_c的计算法如下：

③可见性剔除：

视点较近时的剔除方法：

视点较远时的剔除方法：

f(x，y，z)＝m(x-x_i)+n(y-y_i)+p(z-z_i)＝0 (3.9)

mx+ny+pz＝0 (3.10)

<mrow> <msup> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>&prime;</mo> </msup> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <mi>m</mi> <mrow> <mo>(</mo> <msub> <mi>mx</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>ny</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>pz</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msup> <mi>m</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>n</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>p</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>+</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> </mrow>

<mrow> <msup> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>&prime;</mo> </msup> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <mi>n</mi> <mrow> <mo>(</mo> <msub> <mi>mx</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>ny</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>pz</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msup> <mi>m</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>n</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>p</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>+</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> </mrow>

<mrow> <msup> <msub> <mi>z</mi> <mi>i</mi> </msub> <mo>&prime;</mo> </msup> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <mi>p</mi> <mrow> <mo>(</mo> <msub> <mi>mx</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>ny</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>pz</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msup> <mi>m</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>n</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>p</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>+</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.11</mn> <mo>)</mo> </mrow> </mrow>

计算出所有8个顶点的P′_i′后，取8个坐标点中的极大极小值minX，maxX，minY，maxY，minZ，maxZ；令

x＝maxX-minX

y＝maxY-minY

z＝maxZ-minZ (3.12)

距离系数a的值为：

④缝隙消除与几何过渡：

设每层LOD的起始距离为start_i，i表示当前LOD所在的层级，其所覆盖的距离范围为range，而视点到当前LOD任意一个顶点的距离为dist，则当前的过渡区域[mStart，mEnd]以及该顶点的过渡系数k为：

<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>m</mi> <mi>S</mi> <mi>t</mi> <mi>a</mi> <mi>r</mi> <mi>t</mi> <mo>=</mo> <mrow> <mo>(</mo> <mrow> <msub> <mi>start</mi> <mi>i</mi> </msub> <mo>+</mo> <mi>r</mi> <mi>a</mi> <mi>n</mi> <mi>g</mi> <mi>e</mi> </mrow> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mfrac> <mn>2</mn> <mn>3</mn> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>m</mi> <mi>E</mi> <mi>n</mi> <mi>d</mi> <mo>=</mo> <mrow> <mo>(</mo> <mrow> <msub> <mi>start</mi> <mi>i</mi> </msub> <mo>+</mo> <mi>r</mi> <mi>a</mi> <mi>n</mi> <mi>g</mi> <mi>e</mi> </mrow> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>k</mi> <mo>=</mo> <mn>1.0</mn> <mo>-</mo> <mi>c</mi> <mi>l</mi> <mi>a</mi> <mi>m</mi> <mi>p</mi> <mrow> <mo>(</mo> <mrow> <mfrac> <mrow> <mi>m</mi> <mi>E</mi> <mi>n</mi> <mi>d</mi> <mo>-</mo> <mi>d</mi> <mi>i</mi> <mi>s</mi> <mi>t</mi> </mrow> <mrow> <mi>m</mi> <mi>E</mi> <mi>n</mi> <mi>d</mi> <mo>-</mo> <mi>m</mi> <mi>S</mi> <mi>t</mi> <mi>a</mi> <mi>r</mi> <mi>t</mi> </mrow> </mfrac> <mo>,</mo> <mn>0.0</mn> <mo>,</mo> <mn>1.0</mn> </mrow> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.14</mn> <mo>)</mo> </mrow> </mrow>

其中，clamp(x，a，b)为钳位函数，即

获得过渡系数后，定义顶点所在网格尺寸为n×n，所在的四叉树深度为l，顶点在网格内归一化的坐标为(x，y)，则网格尺寸相对于世界坐标下的实际距离放缩值scale以及过渡后的坐标(x＇，y＇)由如下公式求得：

<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>s</mi> <mi>c</mi> <mi>a</mi> <mi>l</mi> <mi>e</mi> <mo>=</mo> <mfrac> <mn>360</mn> <msup> <mn>2</mn> <mrow> <mi>l</mi> <mo>+</mo> <mn>2</mn> </mrow> </msup> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msup> <mi>x</mi> <mo>&prime;</mo> </msup> <mo>=</mo> <mi>x</mi> <mo>-</mo> <mi>f</mi> <mi>r</mi> <mi>a</mi> <mi>c</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>&CenterDot;</mo> <mfrac> <mi>n</mi> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mfrac> <mn>2</mn> <mi>n</mi> </mfrac> <mo>&CenterDot;</mo> <mi>s</mi> <mi>c</mi> <mi>a</mi> <mi>l</mi> <mi>e</mi> <mo>&CenterDot;</mo> <mi>k</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msup> <mi>y</mi> <mo>&prime;</mo> </msup> <mo>=</mo> <mi>y</mi> <mo>-</mo> <mi>f</mi> <mi>r</mi> <mi>a</mi> <mi>c</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>&CenterDot;</mo> <mfrac> <mi>n</mi> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mfrac> <mn>2</mn> <mi>n</mi> </mfrac> <mo>&CenterDot;</mo> <mi>s</mi> <mi>c</mi> <mi>a</mi> <mi>l</mi> <mi>e</mi> <mo>&CenterDot;</mo> <mi>k</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.15</mn> <mo>)</mo> </mrow> </mrow>

其中，frac(x)表示取x的小数部分；

动态节点的生成：

正则点细分过程：对四边形V₅V₆V₉V₁₀进行细分时，会在其四条边上形成四个边点E₀E₁E₂E₃，同时还有一个面点F，使原来的面片分裂成四个子面片；E₀，E₁，E₂，E₃的值如下公式：

<mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>E</mi> <mn>0</mn> </msub> </mtd> <mtd> <msub> <mi>E</mi> <mn>1</mn> </msub> </mtd> <mtd> <msub> <mi>E</mi> <mn>2</mn> </msub> </mtd> <mtd> <msub> <mi>E</mi> <mn>3</mn> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mi>&alpha;</mi> </mtd> <mtd> <mi>&beta;</mi> </mtd> <mtd> <mi>&beta;</mi> </mtd> <mtd> <mi>&alpha;</mi> </mtd> </mtr> </mtable> </mfenced> <mo>&times;</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>V</mi> <mn>4</mn> </msub> </mtd> <mtd> <msub> <mi>V</mi> <mn>2</mn> </msub> </mtd> <mtd> <msub> <mi>V</mi> <mn>8</mn> </msub> </mtd> <mtd> <msub> <mi>V</mi> <mn>1</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>V</mi> <mn>5</mn> </msub> </mtd> <mtd> <msub> <mi>V</mi> <mn>6</mn> </msub> </mtd> <mtd> <msub> <mi>V</mi> <mn>9</mn> </msub> </mtd> <mtd> <msub> <mi>V</mi> <mn>5</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>V</mi> <mn>6</mn> </msub> </mtd> <mtd> <msub> <mi>V</mi> <mn>10</mn> </msub> </mtd> <mtd> <msub> <mi>V</mi> <mn>10</mn> </msub> </mtd> <mtd> <msub> <mi>V</mi> <mn>9</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>V</mi> <mn>7</mn> </msub> </mtd> <mtd> <msub> <mi>V</mi> <mn>14</mn> </msub> </mtd> <mtd> <msub> <mi>V</mi> <mn>11</mn> </msub> </mtd> <mtd> <msub> <mi>V</mi> <mn>13</mn> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4.1</mn> <mo>)</mo> </mrow> </mrow>

面点F的值由公式4.2得出：

<mrow> <mo>&lsqb;</mo> <mi>F</mi> <mo>&rsqb;</mo> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>E</mi> <mn>0</mn> </msub> </mtd> <mtd> <msub> <mi>E</mi> <mn>1</mn> </msub> </mtd> <mtd> <msub> <mi>E</mi> <mn>2</mn> </msub> </mtd> <mtd> <msub> <mi>E</mi> <mn>3</mn> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>&times;</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mi>&alpha;</mi> </mtd> </mtr> <mtr> <mtd> <mi>&beta;</mi> </mtd> </mtr> <mtr> <mtd> <mi>&beta;</mi> </mtd> </mtr> <mtr> <mtd> <mi>&alpha;</mi> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4.2</mn> <mo>)</mo> </mrow> </mrow>

其中，

②随机分形：用变量(x，y)表示曲面上点的坐标，X(x，y)表示该点在曲面上的高；则FBM曲面定义为在某概率空间上，指数为H的一个随机过程，0<H<1

随机变量X：R→R²满足：

<mrow> <mtable> <mtr> <mtd> <mrow> <mi>P</mi> <mrow> <mo>(</mo> <mi>X</mi> <mo>(</mo> <mrow> <mi>x</mi> <mo>+</mo> <mi>h</mi> <mo>,</mo> <mi>y</mi> <mo>+</mo> <mi>k</mi> </mrow> <mo>)</mo> <mo>-</mo> <mi>X</mi> <mo>(</mo> <mrow> <mi>x</mi> <mo>,</mo> <mi>y</mi> </mrow> <mo>)</mo> <mo><</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>=</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mn>1</mn> <msqrt> <mrow> <mn>2</mn> <mi>&pi;</mi> <msup> <mrow> <mo>(</mo> <mrow> <msup> <mi>h</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>k</mi> <mn>2</mn> </msup> </mrow> <mo>)</mo> </mrow> <mi>H</mi> </msup> </mrow> </msqrt> </mfrac> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <mi>exp</mi> <mrow> <mo>&lsqb;</mo> <mrow> <mo>-</mo> <mfrac> <msup> <mi>r</mi> <mn>2</mn> </msup> <mrow> <mn>2</mn> <msup> <mrow> <mo>(</mo> <mrow> <msup> <mi>h</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>k</mi> <mn>2</mn> </msup> </mrow> <mo>)</mo> </mrow> <mi>H</mi> </msup> </mrow> </mfrac> </mrow> <mo>&rsqb;</mo> </mrow> <mi>d</mi> <mi>r</mi> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4.3</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

Perlin噪声：对于一幅噪声图像，首先给整数坐标点赋予随机方向的法向量与灰度值，除此以外任意一点的像素值由该点所在最邻近的单元格内四个顶点的法向量，和连接该点与四个顶点的方向向量来决定，

设为图中所求点的向量，n为噪声的叠加次数，取n∈[3，12]；设噪声函数波长为λ，则频率ω＝λ^-0.5i，该点的灰度值为：

<mrow> <mi>N</mi> <mi>S</mi> <mrow> <mo>(</mo> <mover> <mi>p</mi> <mo>&RightArrow;</mo> </mover> <mo>)</mo> </mrow> <mo>=</mo> <mi>N</mi> <mrow> <mo>(</mo> <mover> <mi>p</mi> <mo>&RightArrow;</mo> </mover> <msup> <mi>&lambda;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>-</mo> <mi>o</mi> <mi>f</mi> <mi>f</mi> <mi>s</mi> <mi>e</mi> <mi>t</mi> <mo>+</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <mo>&lsqb;</mo> <mi>N</mi> <mrow> <mo>(</mo> <mover> <mi>p</mi> <mo>&RightArrow;</mo> </mover> <msup> <mi>&lambda;</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>-</mo> <mi>o</mi> <mi>f</mi> <mi>f</mi> <mi>s</mi> <mi>e</mi> <mi>t</mi> <mo>&rsqb;</mo> <msub> <mi>E</mi> <mi>i</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4.4</mn> <mo>)</mo> </mrow> </mrow>

其中，offset为偏移量取-0.1，是一个噪声函数；