CN113030898A - 天气雷达基数据立体曲面三角网构建方法 - Google Patents

Abstract

本发明公开了一种天气雷达基数据立体曲面三角网构建方法，S1，数据解析：根据天气雷达基数据存储特征和编码格式对所述基数据进行解析，获取以极坐标表示的三维锥面数据；S2，数据转换：根据指定区域的坐标范围，通过数据插值方法，把以极坐标表示的三维锥面数据转换为笛卡尔坐标表示的三维立体格点场数据；S3，曲面追踪：使用移动立方体算法追踪表示某个要素值的等值曲面，等值曲面由若干个三角形按照一定的拓扑结构连结形成；S4，数据输出：输出保存表示等值曲面的三角网信息。本发明所构建的三维格点场数据提高了垂直方向数据密度，且能够根据指定的雷达探测范围，动态构建三维格点场，使风暴单体或局部范围的三维重构成为现实。

Description

天气雷达基数据立体曲面三角网构建方法

技术领域

本发明涉及天气预报中的探测资料，尤其是涉及天气雷达基数据立体曲面三角网构建方法。

背景技术

新一代多普勒天气雷达探测数据是天气预报中最常用的探测资料之一，对研究降水的形成，分析中小尺度天气系统，警戒强对流天气等具有重要意义，是中小尺度天气监测的重要工具，为气象工作者作出准确的预报提供了很大的帮助。现有对天气雷达回波的分析方法主要停留在二维平面， 比如PPI (平面位置显示) 、CAPPI(等高平面位置显示) 和RHI(距离高度显示) 等平面图像，不能很好的表达风暴云团的三维结构特征。近年来，有学者提出了基于WebGL技术来实现多普勒天气雷达的三维可视化的方法，先对多普勒天气雷达采集生成的基数据文件进行数据预处理，然后采用梯形柱体剖分方法进行等值面的构建，最后利用WebGL技术实现多普勒天气雷达回波的三维图形绘制。如中国专利公开号CN109613540 A、CN105701859A，就公开了利用WebGL技术实现多普勒天气雷达回波的三维图形绘制方法，其主要步骤为：1，读取多普勒天气雷达基数据文件，解析报文获取雷达回波反射率值；2，对获取的雷达回波反射率数据进行坐标转换和灰度数据预处理；3，采用梯形柱体进行三维雷达数据场的等值面构建：选取相邻仰角、相邻方位角和相邻径向距离的8个顶点构成的弧形柱体作为基本处理单元，通过判断顶点回波值和等值面的阈值关系，实现天气雷达三维数据场的不同阈值等值面的逼近表示；4，基于WebGL技术的三维图像渲染和显示,包括三角面片的绘制、颜色和光源的设置、实现透明显示和进行坐标变换实现三维图像的旋转、放大和缩小。

上述使用雷达的多层锥面扫描数据构建三角形面片，采用梯形柱体对整个雷达数据空间进行剖分，存在垂直方向精度不够的不足，导致三维可视化效果不理想。

发明内容

本发明目的在于提供一种基于移动立方体算法的天气雷达基数据立体曲面三角网构建方法，实现快速构建高分辨率的三维等值曲面。

为实现上述目的，本发明采取下述技术方案：

本发明所述天气雷达基数据立体曲面三角网构建方法，包括下述步骤：

S1，数据解析：根据天气雷达基数据存储特征和编码格式对所述基数据进行解析，获取以极坐标表示的三维锥面数据；

S2，数据转换：根据指定区域的坐标范围，通过数据插值方法，把以极坐标表示的所述三维锥面数据转换为笛卡尔坐标表示的规则化的三维立体格点场数据；

S3，曲面追踪：使用移动立方体算法追踪表示某个要素值的等值曲面，所述等值曲面由若干个三角形按照一定的拓扑结构连结形成；

S4，数据输出：输出保存表示等值曲面的三角网信息。

S1中，根据天气雷达基数据存储特征和编码格式，每次读取一条径向数据记录，然后解析分离本条径向数据记录的扫描开始/结束标示符、采样时间、仰角、方位角、距离库数和雷达回波数值，最终得到以距离、方位角和仰角表示的三维锥面数据。

S2中，所述数据转换步骤为：

S2.1, 所述规则化的三维立体格点场设计：

依据雷达探测有效距离和高度，在水平方向上，每公里取一个格点，根据指定的探测范围，动态设置X方向和Y方向格点数；垂直方向取18公里，通过计算机在所述垂直方向上动态设置切片层数，默认36层，每层间隔500米，形成规则化的三维格点场，以提高垂直方向的数据精细化程度；

S2.2, 数据插值：

使用三线性插值把雷达基数据转换为规则化的三维格点场数据；三次线性内插顺序为：先径向内插；再方位内插；最后进行垂直内插，即仰角内插。

S3中，所述曲面追踪步骤为：

S3.1，所述规范化的三维格点场中每个小的立方体单元作为一个体素，所述体素的每个顶点都存在对应的标量值；如果体素顶点上的值大于或等于所述等值面值，则定义该顶点位于等值面之外，标记为“0”； 如果体素顶点上的值小于等值面值，则定义该顶点位于等值面之内，标记为“1”；

S3.2，体素顶点、体素边和三角形面片的定义：

用D₀～D₇标记体素的顶点，用数字0～11标记体素的边；当体素下方的某个顶点的值小于等值面值，其他顶点上的值都大于或等于等值面值，则生成一个与体素的边相交的三角形面片，而所述三角形面片顶点的具体位置则根据等值面值和相应边对应的顶点的要素值，用线性插值计算得到；对于与等值面存在交点的体素边，交点坐标用P表示，P1、P2代表边上两个端点的坐标，v1、v2代表这两个端点上的值，v代表等值面值，则交点坐标P的计算公式如下：

据此计算所有情况下的交点坐标；

S3.3，通过以下三个步骤获得三角形面片数据：

S3.3.1，计算体素顶点状态索引值：

假设V₀～V₇是体素顶点D₀～D₇的要素值，F₀～F₇是所述体素顶点D₀～D₇对应的8位二进制状态标志，定义体素顶点的状态索引cubeindex = F₇F₆F₅F₄F₃F₂F₁F₀，依次检查所述V₀～V₇值，如果V_i值小于等值面值，则F_i置1否则为0，最后得到体素顶点的状态索引值，i=0～7；

S3.3.2，通过顶点状态索引查找并计算三角形顶点坐标；

S3.3.3，计算三角形顶点法向量：

假设体素顶点D(i, j, k)上的值为f(x,y,z)，采用中心差分方法计算该点处的梯度矢量G_x，G_y，G_z：

对所述G_x、G_y、G_z进行归一化，得到顶点D(i, j, k)上的单位法向量，然后使用线性插值得到三角形各个顶点的法向量。

本发明首先对多普勒天气雷达基数据文件进行数据预处理，依据风暴云团的发展高度和雷达回波数据呈锥体分布特点，在垂直方向取18公里，通过计算机动态划分垂直方向切片层，通过三线性插值把雷达基数据转换成规则立方体状的三维格点场数据，然后采用移动立方体算法进行三角网构建，从而生成能够在3D环境中使用的三角形面片数据。本发明所构建的规则的三维格点场数据提高了垂直方向的数据密度，且能够根据指定的雷达探测范围，动态构建三维格点场，使风暴单体或局部范围的三维重构成为现实，弥补了现有方法存在的不足。

附图说明

图1是本发明的流程图。

图2是本发明所述雷达基数据记录格式示意图。

图3是本发明所述规则化的三维立体格点场示意图。

图4是本发明所述三线性插值示意图。

图5是本发明所述移动立方体算法的15种基本情形示意图。

图6a、图6b是本发明所述体素顶点、体素边和三角形面片的定义示意图。

图7是本发明所述交点坐标计算示意图。

具体实施方式

下面结合附图对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述实施例。

如图1所示，本发明所述的天气雷达基数据立体曲面三角网构建方法，包括下述步骤：

S1，数据解析：根据天气雷达基数据存储特征和编码格式对所述基数据进行解析，获取以极坐标表示的三维锥面数据；如图2所示，图2中，N为正北方向，E为正东方向，Z为垂直方向；r为回波点距雷达站距离，θ为雷达扫描仰角，φ为径向波束方位角；h为回波点垂直高度（由r和θ计算得到）。

多普勒天气雷达在多个仰角上有360度方位角的径向探测数据，包含一组相对于雷达站点由近到远不同距离上多个回波点的探测信息，以极坐标方式存储，包括方位角、仰角、采样时间和雷达回波数值。

根据天气雷达基数据存储特征和编码格式，每次读取一条径向数据记录，然后解析分离本条径向数据记录的扫描开始/结束标示符、采样时间、仰角、方位角、距离库数和雷达回波数值，最终得到以距离、方位角和仰角表示的三维锥面数据。

S2，数据转换：根据指定区域的坐标范围，通过数据插值方法，把以极坐标表示的所述三维锥面数据转换为笛卡尔坐标表示的规则化的三维立体格点场数据；步骤为：

S2.1，规则化的三维立体格点场设计：

依据雷达探测有效距离和高度，在水平方向上，每公里取一个格点，根据指定的探测范围，动态设置X和Y方向格点数；垂直方向取18公里，通过计算机在垂直方向上动态设置切片层数，默认36层，每层间隔500米，形成规则化的三维格点场，如图3所示，以提高垂直方向的数据精细化程度。

图3中，三维格点场中任意格点D，其笛卡尔坐标(x,y,z)与极坐标下的距离r、方位角φ、仰角θ的转换公式如下：

（1）

根据公式（1）计算出坐标(x,y,z)对应的极坐标参数(距离r, 方位角φ, 仰角θ)，在极坐标数据场中通过三线性插值即可求得点格点D所在的回波值。

S2.1.2，数据插值：

如图4所示，使用三线性插值把雷达基数据转换为三维格点场数据；三次线性内插顺序为：先径向内插，再方位内插，最后进行垂直内插（仰角内插）。

图4中，f是（r,φ,θ）的函数；设点f(r,φ,θ)落在由f₁(r₁,φ₁,θ₁)、f₂(r₂,φ₁,θ₁)、f₃(r₁,φ₂,θ₁)、f₄(r₂,φ₂,θ₁)、f₅(r₁,φ₁,θ₂)、f₆(r₂,φ₁,θ₂)、f₇(r₁,φ₂,θ₂) 、f₈(r₂,φ₂,θ₂)所围成的锥体内；其中，r为距离，φ为方位角，θ为仰角；r₁是原点o到f₁（f₃、f₅、f₇）的距离,r₂是原点o到f₂(f₄、f₆、f₈)的距离；φ₁是f₁f₂（f₅f₆）所在径向的方位角，φ₂是f₃f₄（f₇f₈）所在径向的方位角；θ₁是面f₁f₂f₄f₃所在扫描锥面的仰角，θ₂是面f₅f₆f₈f₇所在扫描锥面的仰角。设点f₁～f₈的值分别为f⁰ ₁、f⁰ ₂、…f⁰ ₈，目标点f(r,φ,θ)的值为f^v(r,φ,θ),则f^v(r,φ,θ)可由该8个点通过三次线性内插得到，其公式如下：

（2）

公式（2）中，

、

为径向插值权重，

为方位角插值权重，

为仰角插值权重，各权重计算方法如下：

（3）

（4）

（5）

（6）

（7）

（8）

S3，曲面追踪：使用移动立方体算法追踪表示某个要素值的等值曲面，所述等值曲面由若干个三角形按照一定的拓扑结构连结形成，步骤如下：

S3.1，移动立方体算法是一种经典的高分辨率三维曲面构造算法，主要应用于医学领域的可视化场景，例如CT扫描和MRI扫描的3D重建等。本方案中，将该方法用于新一代多普勒天气雷达基数据的立体曲面三角网构建中。

三维格点场中每个小的立方体单元作为一个体素，体素的每个顶点都存在对应的标量值；如果体素顶点上的值大于或等于等值面值，则定义该顶点位于等值面之外，标记为“0”； 如果体素顶点上的值小于等值面值，则定义该顶点位于等值面之内，标记为“1”。

每个体素单元有8个顶点，共存在2⁸ = 256种情形；如图5所示，是移动立方体算法的15种基本情形，其他241种情形可以通过这15种基本情形的旋转、映射等方式实现。

S3.2，体素顶点、体素边和三角形面片的定义：

用D₀～D₇标记体素的顶点，用数字0～11标记体素的边；当体素下方的某个顶点的值小于等值面值，其他顶点上的值都大于或等于等值面值，则生成一个

如图6 a所示，用D₀～D₇标记体素的顶点，用数字0～11标记体素的边；假如体素下方的顶点D₃的值小于等值面值，其他顶点上的值都大于或等于等值面值，则生成一个与体素的边2，3，11相交的三角形面片，如图6b所示；而三角形面片顶点的具体位置则需要根据等值面值和相应边对应的顶点D₃-D₂，D₃-D₀，D₃-D₇的要素值，用线性插值计算得到。

如图7所示，对于与等值面存在交点的体素边，交点坐标用P表示，P1、P2代表边上两个端点的坐标，v1、v2代表这两个端点上的值，v代表等值面值，那么交点坐标的计算公式如下：

（9）

据此即可计算所有情况下的交点坐标。

S3.3，通过以下三个步骤获得三角形面片数据：

S3.3.1，假设V₀～V₇是体素顶点D₀～D₇的要素值，F₀～F₇是顶点D₀～D₇的对应的8位二进制状态标志，定义体素顶点的状态索引cubeindex = F₇F₆F₅F₄F₃F₂F₁F₀，依次检查V₀～V₇，如果V_i（i=0～7）小于等值面值，则F_i置1否则为0，最后得到体素顶点的状态索引值。

S3.3.2，通过顶点状态索引查找并计算三角形顶点坐标：

以图6a、图6b所示为例，仅顶点D₃标记为“1”，其他顶点标记为“0”，那么体素单元的顶点状态索引为8，通过查找表得到TriTable[8] = {3, 11, 2, -1, -1, -1, -1, -1,-1, -1, -1, -1, -1, -1, -1, -1}，其中-1表示无效值；有效的数值（3，11，2）表示三角形面片的3个顶点为体素的边3，11，2和等值面相交的交点；利用公式9即可计算得到三角形的顶点坐标。

S3.3.3，计算三角形顶点法向量：

假设体素顶点D(i, j, k)上的值为f(x,y,z)，采用中心差分方法（公式10）即可计算该点处的梯度矢量G_x，G_y，G_z：

（10）

对梯度矢量G_x，G_y，G_z进行归一化，得到顶点D(i, j, k)上的单位法向量，然后使用线性插值得到三角形各个顶点的法向量。

S4，数据输出：

经过上述3个步骤，生成了每个三角形面片的顶点坐标数据（x,y,z）和顶点单位法向量数据（nx，ny，nz），然后搜索合并相同的顶点坐标并记录其索引号，生成三角形面片的顶点索引数据（idx1,idx2,idx3）。

生成的三角形面片数据中，顶点数据和顶点法向量数据是[-1,1]区间的规范化数据，顶点索引数据是无符号的整型数据；该数据可被应用于不同的3D环境中，比如Direct3D、OpenGL、WebGL，还可以被各种三维建模软件使用。

本发明的术语解释：

天气雷达：我国部署的多普勒天气雷达，是一种主动遥感的探测工具，在测量云、降雨和各种强对流天气发生发展内在因素方面有重要的应用。

基数据：多普勒天气雷达在多个仰角上扫描，每个仰角取360个左右的方位角，每个方位角上得到一条径向回波数据，所有径向回波数据构成一个扫描锥面，称为“一个仰角层”，多个仰角的锥面数据组成一个基数据。在VCP21扫描模式下，每6分钟取得一个基数据，由9个仰角层的锥面数据构成，每层锥面数据包含以雷达站点为原点的多个方位角的径向数据，每条径向数据包含若干个距离库（方位角数量和距离库数量需要根据不同的雷达型号确定），每个距离库表示相对与原点指定距离的回波值。

移动立方体算法：是3D应用领域面绘制的一种经典算法，能够快速构建高分辨率的三维等值曲面。

Claims (4)

1.一种天气雷达基数据立体曲面三角网构建方法，其特征在于：包括下述步骤：

S1，数据解析：根据天气雷达基数据存储特征和编码格式对所述基数据进行解析，获取以极坐标表示的三维锥面数据；

S2，数据转换：根据指定区域的坐标范围，通过数据插值方法，把以极坐标表示的所述三维锥面数据转换为笛卡尔坐标表示的规则化的三维立体格点场数据；

S3，曲面追踪：使用移动立方体算法追踪表示某个要素值的等值曲面，所述等值曲面由若干个三角形按照一定的拓扑结构连结形成；

S4，数据输出：输出保存表示等值曲面的三角网信息。

2.根据权利要求1所述天气雷达基数据立体曲面三角网构建方法，其特征在于：S1中，根据天气雷达基数据存储特征和编码格式，每次读取一条径向数据记录，然后解析分离本条径向数据记录的扫描开始/结束标示符、采样时间、仰角、方位角、距离库数和雷达回波数值，最终得到以距离、方位角和仰角表示的三维锥面数据。

3.根据权利要求1所述天气雷达基数据立体曲面三角网构建方法，其特征在于：S2中，所述数据转换步骤为：

S2.1, 所述规则化的三维立体格点场设计：

依据雷达探测有效距离和高度，在水平方向上，每公里取一个格点，根据指定的探测范围，动态设置X方向和Y方向格点数；垂直方向取18公里，通过计算机在所述垂直方向上动态设置切片层数，默认36层，每层间隔500米，形成规则化的三维格点场，以提高垂直方向的数据精细化程度；

S2.2, 数据插值：

使用三线性插值把雷达基数据转换为规则化的三维格点场数据；三次线性内插顺序为：先径向内插；再方位内插；最后进行垂直内插，即仰角内插。

4.根据权利要求1所述天气雷达基数据立体曲面三角网构建方法，其特征在于：S3中，所述曲面追踪步骤为：

S3.1，所述规范化的三维格点场中每个小的立方体单元作为一个体素，所述体素的每个顶点都存在对应的标量值；如果体素顶点上的值大于或等于所述等值面值，则定义该顶点位于等值面之外，标记为“0”； 如果体素顶点上的值小于等值面值，则定义该顶点位于等值面之内，标记为“1”；

S3.2，体素顶点、体素边和三角形面片的定义：

用D₀～D₇标记体素的顶点，用数字0～11标记体素的边；当体素下方的某个顶点的值小于等值面值，其他顶点上的值都大于或等于等值面值，则生成一个与体素的边相交的三角形面片，而所述三角形面片顶点的具体位置则根据等值面值和相应边对应的顶点的要素值，用线性插值计算得到；对于与等值面存在交点的体素边，交点坐标用P表示，P1、P2代表边上两个端点的坐标，v1、v2代表这两个端点上的值，v代表等值面值，则交点坐标P的计算公式如下：

据此计算所有情况下的交点坐标；

S3.3，通过以下三个步骤获得三角形面片数据：

S3.3.1，计算体素顶点状态索引值：

假设V₀～V₇是体素顶点D₀～D₇的要素值，F₀～F₇是所述体素顶点D₀～D₇对应的8位二进制状态标志，定义体素顶点的状态索引cubeindex = F₇F₆F₅F₄F₃F₂F₁F₀，依次检查所述V₀～V₇值，如果V_i值小于等值面值，则F_i置1否则为0，最后得到体素顶点的状态索引值，i=0～7；

S3.3.2，通过顶点状态索引查找并计算三角形顶点坐标；

S3.3.3，计算三角形顶点法向量：

假设体素顶点D(i, j, k)上的值为f(x,y,z)，采用中心差分方法计算该点处的梯度矢量G_x，G_y，G_z：

对所述G_x、G_y、G_z进行归一化，得到顶点D(i, j, k)上的单位法向量，然后使用线性插值得到三角形各个顶点的法向量。

Images