一种基于狄隆涅三角剖分的雷达反射率数据融合方法

摘要

本发明提供一种基于狄隆涅三角剖分的雷达反射率数据融合方法，包括以下步骤：读取多部雷达的反射率数据信息，进行预处理，提高参与融合反射率数据的质量，形成反射率数据点集；采用虚拟格网划分技术对反射率数据点集进行三维格网划分得到多个四面体；在每一个格网上构建对应的初始四面体，对初始四面体内的反射率数据点采用逐点插入法构建三维delaunay三角剖分；采用空间索引法，遍历各个格网中包含的待插值点，以定位插值点所在的四面体；采用四面体重心坐标插值，并行计算对所有网格点的插值。该方法减少了对同一个待插值点的插值处理次数、避免了多个网格数据融合过程中的叠加，从而减少了误差的累积，提高反射率融合的精度和效率。

权利要求

1.一种基于狄隆涅三角剖分的雷达反射率数据融合方法，其特征在于，包括以下步骤：读取多部雷达的反射率数据信息，进行预处理，提高参与融合反射率数据的质量，形成反射率数据点集；

采用虚拟格网划分技术对反射率数据点集进行三维格网划分；

在每一个格网上构建对应的初始四面体，对初始四面体内的反射率数据点采用逐点插入法构建三维delaunay三角剖分得到多个四面体；

采用空间索引法，遍历各个格网中包含的待插值点，以定位插值点所在的四面体；

采用四面体重心坐标插值，并行计算所有网格点的插值；

其中，所述采用空间索引法，遍历各个格网中包含的待插值点，以定位插值点所在的四面体的步骤包括：

通过一个四面体的四个顶点坐标确定该四面体在Z方向上所占的区间；

判断待插值点所在的水平平面是否位于该区间中，如在，则计算所述水平平面与该四面体的交点；

判断交点的个数，如交点为1，则定位完成；反之，则按照交点的坐标Y值从大到小进行排序，多个交点构成多边形；

计算该待插值点所在水平平面上所有网格点是否位于所述多边形当中：若是，则表明当前四面体包含该待插值点；

否则，开始下一个待插值点的判断；其中，所述插值通过以下方式得到：

获得一个已定位插值点对应的四面体的重心坐标插值；

计算该已定位插值点对应的四面体的各顶点的权重；

根据各顶点的权重以及各顶点的风速计算该已定位插值点的插值；所述插值的计算公式为：

其中，VP为已定位插值点P的插值，VPx、VPy、VPz分别为已定位插值点P在x、y、z方向上的插值，WA、WB、WC、WD分别为已定位插值点P对应所在四面体A、B、C、D四个顶点的权重，分别为顶点A的径向风速在x、y、z方向上的投影， 分别为顶点B的径向风速在x、y、z方向上的投影， 分别为顶点C的径向风速在x、y、z方向上的投影， 分别为顶点D的径向风速在x、y、z方向上的投影，VΔABCD为四面体ABCD的体积，VΔPBCD、VΔPACD、VΔPABD、VΔPABC分别为P与顶点BCD构成的多面体体积、P与顶点ACD构成的多面体体积、P与顶点ABD构成的多面体体积、P与顶点ABC构成的多面体体积。

2.根据权利要求1所述的方法，其特征在于，所述预处理包括对所述反射率数据依次进行平滑滤波、补足缺测点、去除孤立数据点以及去除无效数据点。

3.根据权利要求1所述的方法，其特征在于，将经预处理后的多部雷达的距离库的坐标从以雷达站点为极点的坐标系转换到同一个笛卡尔坐标系下。

4.根据权利要求3所述的方法，其特征在于，所述采用虚拟格网划分技术对预处理后的反射率数据进行三维格网划分的步骤包括：选定目标区域；

将所述目标区域按相同宽度对矩形包围盒进行格网划分，划分成大小近似相等的格网，所述格网为正方形或近似正方形。

5.根据权利要求4所述的方法，其特征在于，所述在每一个格网上构建对应的初始四面体，对初始四面体内的数据采用逐点插入法构建三维delaunay三角剖分得到多个四面体的步骤具体包括：

(1)在每个格网区域全覆盖构建对应的初始四面体，作为每个格网区域对应的delaunay初始四面体；

(2)从待插值的数据点集中选择一个点进行插入操作，并在所有四面体中定位到该点以及该点对应所在的四面体，将该点与所在的四面体的顶点连接，完成点的插入；

(3)依据delaunay三角网的外接空圆特性，利用LOP算法对三维空间中的四面体进行优化；

(4)重复步骤(2)、(3)直到所有反射率数据点集中的点都包含在delaunay三角网中；

(5)将初始四面体和与其连接的所有线条删除。

6.根据权利要求5所述的方法，其特征在于，当交点个数不为1的时候，交点为3或者4个，所述多边形为对应的三角形或四边形。

说明书

一种基于狄隆涅三角剖分的雷达反射率数据融合方法

技术领域

[0001] 本发明属于天气雷达数据处理技术领域，具体涉及一种基于狄隆涅三角剖分的雷达反射率数据融合方法。

背景技术

[0002] 现有的多普勒天气雷达多采用锥形扫描方式工作，这样的方式会使反射率数据密度在空间中分布不均匀，且无法探测到范围内所有点的反射率数据，因此需要将其插值到
感兴趣区域中的均匀的网格点上；而单部天气雷达的扫描范围有限，且存在“静锥区”无法
观测的问题，常利用多部天气雷达联合辅助观测，融合探测资料能够有效弥补这些不足，并
且能够有效提高模式降水预报性能。目前国际上雷达反射率拼图已经成为业务系统数值同
化重要观测资料。

[0003] 若将空间中多个雷达数据看作整体，则失去了单部雷达数据的固定分布特点，因此传统的反射率融合方法常分为两个阶段进行处理，即插值和融合。第一阶段插值是将单
部雷达的反射率资料利用一定的插值方法得到某一高度下的水平网格数据；再将多部雷达
插值得到的多组网格资料进行融合拼接。该方法会将插值过程中带来的误差直接传播到融
合拼接中，相当于对插值后的反射率进行拼图。多个雷达的多普勒数据之间相对独立，没有
直接的联系，误差较大。在各部雷达的插值过程中，常采用垂直线性插值法，垂直水平线性
插值法，自适应Barnes插值法和八点线性插值法等；这些插值方法均是从单个雷达的反射
率数据分布特点为依据，结合待插值网格点相对于雷达站的位置信息，以寻找网格点周围
合适的原始数据点作为影响点。按照线性插值等常用插值方法，将从原始数据中选的一些
影响点进行处理，得到网格点的值。分别将多部雷达探测到的反射率数据插值到同一网格
后，采用距离权重，最大值法或者最邻近法融合多个网格数据。这种方法也可以理解为将雷
达数据从分布不均匀的数据转化为均匀网格数。

[0004] Delaunay(狄隆涅)三角剖分被广泛的应用在数据插值、有限元分析、计算机图形学以及地理信息系统等领域。Delaunay三角剖分就其良好的数学特征，剖分出来的三角形
网格均匀的优化特性而得到了广泛的应用，如有限元分析，计算机图形学等。在二维情况
下，Delaunay三角剖分具有以下特性：唯一性，空圆特性，最大化最小角。通过将散乱的数据
点构成三角网，很好的表达了点与点之间的邻接关系和拓扑结构。值得注意的是，在二维空
间中，散乱的数据经三角剖分后，构成了以三角网为最小单元的网络。在三维空间中，以四
面体为最小的单元网络。这一特点弥补了来自多个观测仪器数据混合后原始数据分布规律
被扰乱的问题。确定待插值点周围的原始数据分布情况，能够便于有效的选取适合的数据
点作为待插值点的影响点。

[0005] 综上所述，现有融合技术的本质是将多个雷达分别插值到同一平面网格后，再进行数据的融合拼接，并未充分利用多部雷达探测使得空间中可参考数据分布密度增加可使
得数据融合精度明显提高的优势；且现有方法分步处理致使工作量繁琐重复，使多次插值
处理和融合过程中误差不断累积。

发明内容

[0006] 有鉴于此，本发明的目的之一在于提供一种基于狄隆涅三角剖分的雷达反射率数据融合方法，该方法能提高反射率融合的精度和效率。

[0007] 为实现上述目的，本发明的技术方案为：

[0008] 一种基于狄隆涅三角剖分的雷达反射率数据融合方法，包括以下步骤：

[0009] 读取多部雷达的反射率数据信息，进行预处理，提高参与融合反射率数据的质量，形成反射率数据点集；

[0010] 采用虚拟格网划分技术对反射率数据点集进行三维格网划分；

[0011] 在每一个格网上构建对应的初始四面体，对初始四面体内的反射率数据点采用逐点插入法构建三维delaunay三角剖分得到多个四面体；

[0012] 采用空间索引法，遍历各个格网中包含的待插值点，以定位插值点所在的四面体；

[0013] 采用四面体重心坐标插值，并行计算对所有网格点的插值。

[0014] 进一步地，所述预处理包括对所述反射率数据依次进行平滑滤波、补足缺测点、去除孤立数据点以及去除无效数据点。

[0015] 进一步地，还包括步骤，将经预处理后的多部雷达的距离库的坐标从以雷达站点为极点的坐标系转换到同一个笛卡尔坐标系下。

[0016] 进一步地，所述采用虚拟格网划分技术对预处理后的反射率数据进行三维格网划分的步骤包括：

[0017] 选定目标区域；

[0018] 将所述目标区域按相同宽度都矩形包围盒进行网格划分，划分成至少大于1个大小相等的格网，所述格网为正方形或近似正方形。

[0019] 进一步地，所述在每一个格网上构建对应的初始四面体，对初始四面体内的数据采用逐点插入法构建三维delaunay三角剖分得到多个四面体的步骤具体包括：

[0020] (1)在每个格网区域内构建对应的初始四面体，作为每个格网区域对应的delaunay初始四面体；

[0021] (2)从待插值的数据点集中选择一个点进行插入操作，并在所有当前的四面体中定位到该点以及该点对应所在的四面体，将该点与所在的四面体的顶点连接，完成点的插
入；

[0022] (3)依据Delaunay三角网的外接空圆特性，利用LOP算法对当前所述四面体进行优化；

[0023] (4)重复步骤(2)、(3)直到所有反射率数据点集中的点都包含在所述四面体中；

[0024] (5)将初始四面体和与其连接的所有线条删除

[0025] 进一步地，所述采用空间索引法，遍历各个格网中包含的待插值点，以定位插值点所在的四面体的步骤包括：

[0026] 通过一个四面体的四个顶点坐标确定该四面体在Z方向上所占的区间；

[0027] 判断待插值点所在的水平平面是否位于该区间中，如在，则计算所述水平平面与该四面体的交点；

[0028] 判断交点的个数，如交点为1，则定位完成；反之，则按照交点的坐标Y值从大到小进行排序，多个交点构成多边形；

[0029] 计算该待插值点所在水平平面上所有网格点是否位于所述多边形当中：

[0030] 若是，则表明当前四面体包含该待插值点；

[0031] 否则，开始下一个待插值点的判断。

[0032] 进一步地，当交个个数不为1的时候，交点为3或者4个，所述多边形为对应的三角形或四边形。

[0033] 进一步地，所述插值通过以下方式得到：

[0034] 获得一个已定位插值点对应的四面体的重心坐标插值；

[0035] 计算该已知定位插值点对应的四面体的各顶点的权重；

[0036] 根据各顶点的权重以及各顶点的风速计算该已定位插值点的插值。

[0037] 进一步地，所述各顶点的权重为已定位插值点所对应的四面体的体积与四面体的体积之间的比值。

[0038] 进一步地，所述插值的计算公式为：

[0039]

[0040] 其中，VP为已定位插值点P的插值，VPx、VPy、VPz分别为已定位插值点P在x、y、z方向上的插值，WA、WB、WC、WD分别为已定位插值点P对应所在四面体A、B、C、D四个顶点的权重，
分别为顶点A的径向风速在x、y、z方向上的投影， 分别为顶点B的
径向风速在x、y、z方向上的投影， 分别为顶点C的径向风速在x、y、z方向上的
投影， 分别为顶点D的径向风速在x、y、z方向上的投影，VΔABCD为四面体ABCD的
体积，VΔPBCD、VΔPACD、VΔPABD、VΔPABC分别为P与顶点BCD构成的多面体体积、P与顶点ACD构成的
多面体体积、P与顶点ABD构成的多面体体积、P与顶点ABC构成的多面体体积。

[0041] 有益效果

[0042] 本发明提供一种基于狄隆涅三角剖分的雷达反射率数据融合方法，该方法本不同于传统两步走的方法，减少了对同一个待插值点的插值处理次数、多个网格数据融合过程
中的叠加，从而减少了误差的累积，提高反射率融合的精度和效率；该方法在插值过程中实
现了并行处理，提高了运行效率；同时该方法不仅能结合多普勒天气雷达观测到的反射率
数据，还能同时结合其他天气观测仪器共同使用，例如相控阵雷达、全固态雷达等，结合多
种雷达观测资料及其他探测资料进行数据融合，从而辅助对天气系统发生发展做出判断
据。

附图说明

[0043] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍。显而易见地，下面描述中的附图是本发
明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以
根据这些附图获得其它的附图。

[0044] 图1为本发明一种基于狄隆涅三角剖分的雷达反射率数据融合方法的一实施例流程图；

[0045] 图2为本发明的狄隆涅三角剖分中数据点分布示意图；

[0046] 图3为本发明的狄隆涅三角剖分效果图；

[0047] 图4为本发明实施例提供的基于四面体重心插值结构示意图；

[0048] 图5为自适应Barnes插值法在3Km高度层上的反射率数据融合效果示意图；

[0049] 图6为本发明中基于狄隆涅三角剖分的四面体重心插值法在3Km高度层上的反射率数据融合效果示意图；

[0050] 图7为自适应Barnes插值法在4Km高度层上的反射率数据融合效果示意图；

[0051] 图8为本发明中基于狄隆涅三角剖分的四面体重心插值法在4Km高度层上的反射率数据融合效果示意图；

[0052] 图9为发明实施例提供的3Km高度处反射率融合数据连续性示意图；

[0053] 图10为本发明实施例提供的在构建的狄隆涅三角剖分网中并行插值得到的斜面数据融合效果示意图。

具体实施方式

[0054] 为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例是
本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员
在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0055] 所举实施例是为了更好地对本发明进行说明，但并不是本发明的内容仅局限于所举实施例。所以熟悉本领域的技术人员根据上述发明内容对实施方案进行非本质的改进和
调整，仍属于本发明的保护范围。

[0056] 实施例1

[0057] 参考图1为一种基于狄隆涅三角剖分的雷达反射率数据融合方法，包括以下步骤：

[0058] S100：读取多部雷达的反射率数据信息，进行预处理，提高参与融合反射率数据的质量，形成反射率数据点集；然后执行步骤S200；

[0059] 本实施例中，接收多部雷达的反射率数据，雷达探测的反射率数据中，有部分数据为无效数据，因此在进行数据融合之前需要删除无效数据，并进行数据预处理，预处理包括
对反射率数据依次进行平滑滤波、补足缺测点、去除孤立数据点以及去除无效数据点，以提
高各部雷达探测到的反射率数据质量；然后将大范围的各部雷达数据看作为一个整体，并
把各距离库(雷达反射率数据在径向上的最小分辨单元，每一个距离库均有其反射率值)的
坐标从以雷达站点为极点的坐标系转换到同一个笛卡尔坐标系下；即将各个雷达的所有距
离库从极坐标系(以各自雷达站点为极点)转换到同一个笛卡尔坐标系下，因为来自不同雷
达的这些距离库所在的坐标系是不一样的，本发明中通过坐标转换，将所有距离库统一到
同一坐标系下。

[0060] S200：采用虚拟格网划分技术对反射率数据点集进行三维格网划分；然后执行步骤S300；

[0061] 雷达探测数据量通常较大，本实施例中采用虚拟格网划分技术对预处理后的反射率数据进行三维格网划分，以便于提高四面体的构建和定位的效率；

[0062] 本步骤中，首先选定目标区域，然后把整个区域按一定的宽度对矩形包围盒进行格网划分，划分成若干大小近似相等的格网，保证每一个格网单元为正方形或近似正方形。

[0063] 优选地，在参与构建delaunay四面体网络的数据点数目较多时，即反射率数据点集中数据量大时，对所有数据进行格网划分，划分时需要保证各个格网区域内的数据应大
于一定的阈值，本发明中的阈值可根据雷达数据的分布密度情况，自行设置。在现有常见的
逐点插入法构建三角网时，对于点定位的算法是遍历三角形链表，找出点所在的三角形；当
存在大量点的加入时，三角形的数量会逐渐增加，效率很低，尤其是在三维空间构建四面体
网时，其效率受到的影响更大。建立虚拟格网将原始数据点集进行分区域处理，不仅提高了
整体数据构建delaunay四面体网络的效率，还有效促进了点定位的处理效率。

[0064] S300：在每一个格网上构建对应的初始四面体，对初始四面体内的反射率数据点采用逐点插入法构建三维delaunay三角剖分得到多个四面体；然后执行步骤S400；

[0065] 本实施例中，首先在步骤S200中划分的每个格网区域全覆盖构建对应的初始四面体，该初始四面体需要足够大到包围该网格内的所有反射率数据点集，作为每个格网区域
对应的delaunay初始四面体；然后完成反射率数据点集中所有数据的插入；

[0066] 具体地，在反射率数据点集中选择一个点，并在所有四面体中定位到该点以及该点对应所在的四面体，将该点与所在四面体的所有顶点连接，完成点的插入；在一般情况
下，Delaunay三角剖分法被使用在二维空间，其构建的三角网为最小单元，而在本发明中，
拓展到三维空间，构建的最小单元即为四面体，因此，本实施例也参考二维空间中Delaunay
三角网的外接空圆特性，利用LOP算法对三维空间中的四面体进行优化，通过逐步优化网
格，减少畸形三角单元的产生，保证高质量的网格生成，该点的数值插入完成；接着重复上
述方法，将反射率数据点集中的所有点都插入构建三维delaunay三角剖分，使得所有反射
率数据点集中的点都包含在Delaunay三角网；最后删除初始四面体和与初始四面体连接的
直线，本实施例中的图2为狄隆涅三角剖分中数据点分布示意图，图3为经过步骤S300 
delaunay四面体剖分得到的四面体网络。

[0067] S400：采用空间索引法，遍历各个格网中包含的待插值点，以定位插值点所在的四面体；然后执行步骤S500；

[0068] 本实施例中格网是针对于已有数据，将其进行划分，划分出来形成虚拟区域称为格网；而网格点即是在待插值面的待插值的点，这些点取自于想要进行数据融合区域上自
己设置的虚拟格点，一般呈均匀分布，以表征目标面的情况；

[0069] 本步骤中，通过空间索引技术定位插值点(雷达反射率数据点集中的点)所属哪一个四面体，实现插值点的定位。常用的定位思想是针对于一个插值点，遍历所有已构建好的
四面体，判断是否包含当前插值点。这种遍历方法增加了运行成本。在实际情况中，待插值
点常以平面上一系列均匀划分并且规则的网格点的形式出现，待插值点集之间具有较好的
联系。因此本发明通过遍历四面体，找到四面体包含哪些插值点。对于一个待插值点只属于
一个四面体，通过依次取出各个四面体，计算所有落在该四面体内的待插值点，遍历完所有
四面体即完成了所有待插值点的定位。

[0070] 具体地，通过一个四面体的四个顶点坐标确定该四面体在Z方向上所占的区间S；然后判断待插值点所在的水平平面是否位于该区间S中，如在，则计算水平平面与该四面体
的交点；随后判断交点的个数，如交点为1，则定位完成；反之，多个交点构成多边形，即当交
点个数不为1的时候，交点为3或者4个，多边形为对应的三角形或四边形则按照交点的Y值
从大到小进行排序；然后计算该待插值点所在水平平面上所有网格点是否位于多边形当
中；对所有的四面体重复遍历以上步骤，直到所有待插值点遍历完。

[0071] S500：采用四面体重心坐标插值，并行计算对所有网格点的插值。

[0072] 本实施例中，假设已定位插值点P处于四面体ABCD的内部，该插值点P与四面体ABCD的任三个点可获得4个多面体，然后获得一个已定位插值点P与这个4四个多面体的重
心坐标插值，可计算该已知定位插值点对应的四面体的各顶点的权重，其基于四面体重心
插值结构示意图可参考图4，本实施例中的各顶点的权重为已定位插值点所对应的四面体
的体积与四面体的体积之间的比值：

[0073]

[0074]

[0075]

[0076]

[0077] WA、WB、WC、WD分别为已定位插值点P对应所在四面体A、B、C、D四个顶点的权重，VΔABCD、VΔPBCD、VΔPACD、VΔPABD、VΔPABC分别为四面体ABCD的体积、P与顶点BCD构成的多面体体积、
P与顶点ACD构成的多面体体积、P与顶点ABD构成的多面体体积、P与顶点ABC构成的多面体
体积。

[0078] 然后根据各顶点的权重以及各顶点的风速计算该已定位插值点的插值：

[0079]

[0080] 其中，VP为已定位插值点P的插值，VPx、VPy、VPz分别为已定位插值点P在x、y、z方向上的插值，WA、WB、WC、WD分别为已定位插值点P对应所在四面体A、B、C、D四个顶点的权重，
分别为顶点A的径向风速在x、y、z方向上的投影， 分别为顶点B的
径向风速在x、y、z方向上的投影， 分别为顶点C的径向风速在x、y、z方向上的
投影， 分别为顶点D的径向风速在x、y、z方向上的投影。

[0081] 实施例2

[0082] 本实施例中，对本发明基于狄隆涅三角剖分的雷达反射率数据融合方法的应用效果作详细说明：

[0083] 本实施例采用现常用的反射率数据插值算法之一：自适应Barnes插值法与本发明的狄隆涅三角剖分算法进行对比实验；参考图5‑图8，其中，图5与图7分别为自适应Barnes
插值法在3Km、4KM高度层上的反射率数据融合效果示意图；图6与图8为分别本发明中基于
狄隆涅三角剖分的四面体重心插值法在3Km、4KM高度层上的反射率数据融合效果示意图，
根据这四幅图可知，其反射率融合结果在这两个高度的效果相差不大，具有基本相同的特
征轮廓；但从总体上来讲，狄隆涅三角剖分算法的融合后数据整体大于传统方法，本实施例
中还对3Km高度处某一直线的融合数据进行分析，其结果可参考图9可以看出，虽然，数据变
化趋势一致，起伏幅度大致相同，基于狄隆涅三角剖分算法的重心插值法能够对雷达反射
率数据进行有效的融合。

[0084] 此外，本发明提出的雷达反射率数据融合方法能够在已构建的三维狄隆涅四面体网络的基础上，多次任意设置待插值剖面上的网格点并快速经过插值得到剖面数据；例如
图10给出了一倾斜剖面经过数据融合后的效果图。

[0085] 上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员
在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多
形式，这些均属于本发明的保护之内。