天气雷达三维回波实时渲染算法研究

王 兴， 钱代丽， 朱 彬，， 安俊琳， 吕晶晶

（南京信息工程大学 a.大气科学与环境气象国家级实验教学示范中心，b.大气科学与气象信息国家级虚拟仿真实验教学中心，南京 210044）

0 引 言

曲面重建是对离散的空间数据进行三维渲染的主要技术手段，其基本思想是按照一组给定的阈值从三维点云数据中提取相应的等值面，利用曲面渲染算法生成一组由不同颜色相互叠加的等值面［5-6］。曲面重建的算法很多，主要可分为两类：一是基于断层轮廓线的表面重建，即将三维数据分割成若干个二维切面，分别计算各切面的等值线，再将这些等值线连接成等值曲面；二是基于体素的等值面重建，即对三维点云构成的体数据场中的各个体素加以分析，提取等值面信息，再将这些等值面连接成等值曲面［7-8］。其中，移动立方体法（Marching Cubes，MC）［9］是沿用至今且颇具影响力的体素等值面重建算法之一。

曲面重建技术在医疗、建筑等工程领域的研究最为深入，近年来，也有一些研究将该技术应用于天气雷达、卫星和数值模式等资料，实现气象数据的三维可视化交互［10-11］。由于雷达探测数据的空间分布极不规则，距离雷达近端与远端的相邻离散点之间，间隔距离相差超过100倍，在重构均匀三维格点场雷达回波时，传统空间插值算法并不适用。雷达远端的插值点不仅难以还原真实的探测值，还极大地增加了点云的密度，造成面片的计算规模陡增，网络结构复杂化，大幅增加了三维渲染的时长，影响三维可视化的实时性和可交互性。

为解决上述技术难题，针对多普勒天气雷达探测方式及其数据特征，在移动立方体算法基础上，提出一系列改进，并对改进前、后的渲染效果做出比对分析。

1 基本原理与技术路线

1.1 移动立方体（MC）算法

MC算法最早是Lorensen等［9］提出，用于解决三维物体的面绘制问题。其本质是将三维坐标系下的物体分解为体素集合，每个体素包含8个空间上紧相邻的顶点，构成一个正立方体，将该正立方体视为基本单位。曲面重建时，先设定一组大小不同的阈值，通过计算每个体素内的梯度变化，并与各个阈值分别进行比较，找出所有被阈值穿过的正立方体，再利用插值算法计算出这些表面（即等值面），再通过绘制一系列阈值的等值面，实现三维矢量图像的绘制。

受制于20世纪80年代有限的计算机图形处理能力，MC算法巧妙地利用三角面片近似表示待求取的等值面，极大地提升了计算效率。由于等值面与体素边界面的交线是双曲线，在特定情况下对相同的等值点可以采取不同的连接方式，这种近似表示的算法无法保证三角片所构成的等值面的拓扑一致性，造成很多三角面片衔接处出现孔隙和不连续现象。为克服这一问题，很多学者对MC算法提出了改进，有的是利用插值和双曲渐近线将体素进一步拆分［12-13］，有的是将MC的正立方体进一步简化，采用四面体剖分［14］等。

1.2 MC应用于雷达回波渲染的技术缺陷

气象雷达探测得到的数据，其空间分布是一个不规则的三维点云结构，随着探测距离的增加，点云愈发稀疏，距离雷达近端与远端的相邻离散点之间，间隔距离相差可达100倍以上。在重构均匀三维各点场雷达回波时，诸如双线性插值、最近邻插值、三次内插等传统空间插值算法并不适用。

宴会中，马老当场宣布，这次活动的成功举办，在座的各位功不可没。公司将论功行赏，给予在座的各位同仁以物质奖励，对公司的高管层将另行奖赏。奖励田卓小姐、高潮先生赴巴厘岛六日游，冯可儿小姐本来也在去巴厘岛之列，但考虑到公司不可一日无主，就留下来代理田卓行使管理权，但公司会额外补偿冯可儿小姐不能成行的损失。

尽管MC利用三角面片拟合曲面的算法设计非常巧妙，在曲面渲染时也表现出不错的性能，但对于面片的生成，却需要消耗“过多”的内存和计算时长。这主要是因为MC算法需要遍历三维点云数据中的所有体素，而与等值面相交的体素数量仅占所有体素的1/10甚至更低。雷达探测的空间范围大（数百公里），数据密度高，一次计算需要处理数百万个体素，而近9成的时间被用于计算这些无效体素。雷达远端的插值点不仅难以还原真实的探测值，还极大地增加了点云的密度，使得体素产生的面片面积过小，造成面片的计算规模陡增，网络结构复杂化，大幅增加了三维渲染的时长，影响了三维回波渲染的实时性和可交互性。

1.3 雷达三维回波实时渲染算法流程

为克服上述不足，本算法以传统MC算法为基础，针对雷达基数据的空间结构特点，着重在三角面片顶点位置的计算、三角剖分构型唯一性的判定以及包含等值面体素的加速遍历，这3方面提出修改和优化。雷达三维回波实时渲染的总体算法流程如图1所示。

图1 总体算法流程图

2 雷达三维回波实时渲染算法

2.1 数据准备

国内常用的天气雷达分为S波段、C波段和X波段3类，每个波段各有多种型号，每个型号的雷达又可采用多种方式的体扫模式。以我国东南沿海地区最常用的SA型为例，它的体扫模式主要有VCP11（Volume Coverage Pattern）、VCP21和VCP31这3种。其中，VCP21常被用于探测非强对流的明显降水情况，该模式的体扫周期为6 min，方位角方向的分辨率为1°，径向分辨率为0.25 km，探测距离为460 km［15］。1个体扫周期内完成9个不同仰角扫描，仰角与高度、水平距离的关系如图2所示。

图2 SA雷达VCP21体扫模式

图2中：hd为水平距离；vh为垂直高度，单位均为km；ev为探测仰角，单位为（°）。通常情况下，MC算法分析的数据在空间上等间隔分布结构，即各个位置的体素密度相同。显然，雷达探测数据结构并不满足这一特征，通行的做法是将极坐标系下的数据结构通过空间插值转换成平面直角坐标系下的三维结构，但额外增加的体素信息不能准确地反映大气中的水气状况，除了方便计算，并不能带来性能上的有益提升。本算法不做插值和坐标系转换，而直接根据雷达探测所构成的数据空间结构，分析各个等值面阈值与基本反射率之间的数值关系。

如图3所示，a～h为雷达探测可识别分辨率下某处相邻的8个顶点，各顶点的数值为该位置探测到的基本反射率，用Z（λ，θ，φ）表示，单位为dBZ。这8个顶点构成一个六面体，将其视作一个体素，它是后续MC计算的基本单位。

图3 极坐标系下相邻8顶点的空间分布

2.2 极坐标系下三角面片顶点位置的计算

体素中的三角面片是拟合等值面的基本单位，要得到三角面片的空间位置，就需要准确计算等值面在体素棱边相交点的坐标。通常在由X、Y、Z构成的空间直角坐标系下，可以利用线性插值计算各个交点的坐标。雷达数据是极坐标系下的空间三维结构，计算方法需做相应转换。假设D1=（λ1，θ1，φ1）和D2=（λ2，θ2，φ2）是体素棱边上两个端点的坐标，阈值为T1的等值面在该棱边上有交点，那么该交点的坐标P（λ，θ，φ）计算方法为：

2.3 基于相邻体素状态的快速判别算法

VCP21体扫模式类似图3所示，由8个顶点构成的六面体（体素）约有1.32×106个，采用MC算法对所有体素逐一遍历，将消耗大量时间，且大多数体素与待求取的等值面并不相交。为提升计算效率，采用一种区域拓扑扩展的思想，先任意选取一个与等值面相交的体素，将其作为基准，筛选出与该体素相邻且包含等值面的体素，再分别将这些体素作为基准，采用相同方法扩展筛选，直到找出所有与等值面相交的体素，再做后续处理。

上述过程的关键是如何判别与基准体素相邻的体素是否与等值面相交。如图4所示，有4个相邻的体素，每个体素包含8个顶点，其中CD是4个体素公有的边，ABCD是上方2个体素公有的面。如果A、B、C、D4个顶点的基本反射率数值不全部小于或大于等值面的阈值，则可判定基准体素的相邻体素与等值面相交。图4中MN表示左上角的基准体素与等值面的交线，由于MN也存在于右上角的体素中，因此，右上角的体素也与等值面相交。

图4 相邻体素等值面与交线的示意图

尽管A、B、C、D的取值不定，但相较于等值面的阈值，只有小于、大于或等于这3种状态。为简化表述，不妨假设当顶点的基本反射率数值小于或等于等值面的阈值时，设定为状态“T”，当大于时设定为状态“F”。那么，对于体素中任意一面的状态可总结为如图5所示的4种状态。

图5 体素中一个面与等值面相交的情况

图5中，体素一个面的顶点用黑色实心圆和白色空心圆分别表示状态“T”和“F”。其中，图5（a）有1个顶点与其他顶点的状态不同；图5（b）、（c）各有2个顶点与其他顶点的状态不同；图5（d）中4个顶点为同一状态。图5中的虚线为体素的一个面与等值面相交的示意，相交的准确位置可采用2.2节所述方法计算得到。其他情况均可通过旋转或翻转图形归并到图5所示的4种情况之一。在遍历计算相邻体素时，可借助哈希表辅助算法执行［16］。

2.4 体素中三角剖分构型的唯一性判定算法

传统MC算法中，当状态“T”和“F”的顶点如图5（c）所示分别位于对角线的两端，那么等值面与体素相交的方式会有两种，因此存在歧义，如图6所示。为保证三角面片拓扑结构的一致性，一些学者提出了子构型拆分和双曲线渐近线判别［17］等方法。无论是子构型拆分中的递归迭代，还是双曲线渐近线交点的三线性插值计算，都使得计算的时间复杂度增加了一个量级。

图6 等值面与体素相交方式的二义性

为加快体素等值面的判定，并消除二义性，本算法借鉴双曲线渐近线判别法的思路［18］，并提出一点改进。通常情况下，等值面与体素相交形成的交线是一对双曲线，该双曲线与体素面上的位置关系如图7所示。

图7 双曲线与体素面上的位置关系

可见，图7（a）中双曲线的两支同时与体素的面相交，双曲线将体素的面分割成3个区域。分别将两组对边上的交点用虚线连接，这两条虚线形成的交点一定和其中一条对角线上的两个顶点同处一个区域，即同为状态“T”或“F”。通过比较该交点处的基本反射率值与阈值的大小关系，即可判定该二义性面与等值面的连接方式。

用图8举例说明，P1～P4是等值面与体素中一个面相交形成的4个交点，这4个交点的坐标可用2.2节所述方法计算得到，再通过联立P1、P2和P3、P4的直线方程，可计算出两条直线相交处O的坐标。由于A～D4个顶点的基本反射率值是已知的，因此可通过双线性插值，计算出O点的基本反射率值，计算方法为：

图8 二义性面的唯一性判定

式中：λA、θA和φA分别表示顶点A极坐标下的3个变量；Z（A）顶点A处的基本反射率值；Z（O）即为两条直接相交处O点的基本反射率值。同样用状态“T”和“F”来表示该值与等值面阈值的大小关系。当Z（O）的状态为“F”时，可确定图8（a）的构型，当当Z（O）的状态为“T”时，可确定图8（b）的构型。

3 实验结果及分析

3.1 实验结果

为实现雷达三维回波的可视化，需将上述过程计算得到的一系列等值面用特定格式的文件输出。目前主流的格式有：COLLADA是一种基于XML的3D交互文件格式，它通过对场景、节点、几何形状、着色器、力学和运动过程的描述，实现对三维矢量图形的渲染和交互。图形语言传输格式（Graphics Language Transmission Format，GLTF）是另一种常用的图形语言交换格式，它用比XML更简洁的JSON格式存储信息，对OpenGL更友好，支持大数据量的3D Web实时渲染。

本实验选用GLTF格式，将一组阈值的等值面依次输出，通过Blender软件加载20 dBZ等值面渲染后的效果如图9所示。

图9 采用Blender渲染的雷达三维回波（20 dBZ）

在天气分析等教学实践过程中，往往需要结合地理信息对雷达回波的特征加以分析，目前诸如Three.js、Cesium和PlayCanvas等开源软件都对GLTF格式数据有很好的支持［19］。其中，Cesium是一个基于JavaScript可创建三维场景和虚拟地球的地图可视化框架，在基于Web的地图动态数据可视化方面具有优异的性能。

本实验将上述阈值的等值面依次加载到Cesium，并用气象上常用的色标对不同阈值的等值面渲染相应颜色，叠加地球影像和地图图像的显示的效果如图10所示。图10（a）、（b）加载的是一同GLTF文件，不同之处在于观测的视角和加载的地图信息。

图10 采用Cesium渲染的雷达三维回波

3.2 结果分析

改进传统MC算法的计算效能，提升雷达三维回波渲染和交互的实时性是本算法改进的主要目标。为得到客观量化的评价指标，以下测试采用相同的软硬件环境。主要硬件配置为：Intel Xeon E5-1630 v4CPU，默认主频3.7 GHz，DDR3 2400 MHz 128 GB内存。算法性能的评价主要包括两个方面：一是计算耗时，它是评估等值面提取效率的重要指标，直接关系到是否满足雷达数据实时渲染的性能要求；二是所生成的三角面片的总数量，虽然更多的三角面片使得图形更加平滑，能有效减少锯齿，但对内存的开销很大。

表1中，“改进1”是以传统MC算法为基础，采用2.1、2.2节所述方法，不做空间插值和坐标转换；“改进2”是在“改进1”的基础上，增加2.3节所述方法进行相邻体素状态的快速判别。针对传统MC的等值面二义性问题，上述测试均采用了2.4节所述方法进行处理。测试过程选用了不同时刻的雷达基数据，分别记录不同算法下的计算耗时和三角面片数量，表1中的数值是10次测试的算法平均值取整。可见，传统MC完成一次渲染需要27 s，其中大量时间消耗在坐标转换过程的插值计算。以及等值面提取时对大量无价值体素的分析。“改进1”优化了上述缺陷，既降低了计算耗时，也减少了三角面片的数量，通过Edge浏览器加载Cesium和GLTF数据，监测内存发现，Edge进程占用内存减少了56%。“改进2”进一步降低了体素遍历的数量，将计算耗时再度提升34%。

表1 改进算法性能对比

在实际操作过程中，采用“改进2”的计算方法，每次加载新的雷达基数据仍需近7 s的等待响应时间，如果将这些雷达基数据提前预处理好，生成GLTF文件存储于服务器中，那么网络加载GLTF文件并通过浏览器渲染的耗时将大幅降低到500 ms以内，达到实时交互、快速响应的目标。通过网页异步传输、预加载等技术，可实现时序三维回波动画的平滑播放，极大增强师生在天气分析教学活动中的体验。

4 结 语

利用曲面重建技术实现天气雷达三维回波的实时渲染，为雷达气象学、基于雷达资料的短临天气分析等教学实践活动提供了更加直观、立体、可交互的应用软件平台。本文设计的算法，针对雷达基数据的空间结构特点，着重在三角面片顶点位置的计算、三角剖分构型唯一性的判定以及包含等值面体素的加速遍历3方面提出修改和优化，测试结果表明，改进后的算法计算耗时缩短至25%，内存占用降低了56%，满足师生教学过程中，雷达数据渲染与交互的实时性要求。该算法也可应用于数值天气模式、气象卫星等资料的三维渲染，具有一定的推广应用价值。