CN102117227A - 天气雷达数据的多核并行计算方法 - Google Patents

Abstract

一种天气雷达数据的多核并行计算方法，其特征是它包括以下步骤：首先，将探测到的天气雷达观测数据用球面坐标系进行存储；其次，将存储的球面坐标系数据进行格点化处理，也就是把球面坐标系下的空间分辨率不均匀的雷达数据插值到统一的笛卡尔坐标下形成空间分辨率均匀的互不相关的网格点数据；第三，利用OPENMP程序将得到的网格点数据在各自线程的控制下进入对应的CPU核中运行即可得到相应的气象数据；最后，将各CPU核计算所得数据进行合成即可得到所需的气象数据。本发明的多核并行技术可使CPU的最大使用率达到97％，充分使用了计算机的运算能力。提高气象数据计算分析的实时性。

Description

天气雷达数据的多核并行计算方法

技术领域

本发明涉及一种气象雷达数据处理方法，尤其是一种利用多核计算机实时处理雷达扫描数据的方法，具体地说是一种天气雷达数据的多核并行计算方法。

背景技术

众所周知，计算机CPU在单核心的的频率达到近乎极限以后，CPU进入多核时代。所谓的多核是指在一枚处理器中集成两个或者多个完整的计算引擎。

多核编程，很重要的一点就是要将多个CPU核同时运转起来，提高CPU的利用率，也就是要提高多核CPU的加速比系数。

但是并不是任何时候都可以同时将多个核运转起来，比如一些不能并行化执行的代码，特别是使用了锁的代码，它的执行是串行化的，无法同时利用多个核，还有负载平衡问题，多个任务无法划分得很均匀，导致实际情况中有些任务先运行完，有些任务后运行完，先运行完的CPU核就处于空闲状态。

多核计算的运用场合是针对大量数据并需要针对大量的数据需要复杂的计算的情况，在数据量很大且对数据的运算量很大的气象领域极为适用。

在中尺度气象领域，由于中尺度现象无论在空间上还是在时间上的分辨率都很高，要求业务数据预报模式或者数值研究模式在空间上和时间上的分辨率尽量高，模式运行的时间尽量短。因此对计算机的计算速度要求很高，多核并行计算是一条很好的解决方案。

新一代天气雷达具有很高的时间和空间分辨率，尤其双线偏振雷达，所能探测到的数据量数倍于常规的多普勒雷达，对超大规模的天气雷达数据的处理，计算的复杂程度比较大，计算的时间开销显著。因此，对于这种大规模数据进行并行处理就显得十分重要。

发明内容

本发明的目的是针对目前的雷达数据处理过程中由于待处理数据的关联性较强，造成各CPU核运行的忙闲不均，时间开销大的问题，发明一种先将雷达数据格点化成互不相关的数据送入各自的CPU核独立运行，最后进行合成得到所需相关气象数据的天气雷达数据的多核并行计算方法。

本发明的技术方案是：

一种天气雷达数据的多核并行计算方法，其特征是它包括以下步骤：

首先，将探测到的天气雷达观测数据用球面坐标系进行存储；

其次，将存储的球面坐标系数据进行格点化处理，也就是把球面坐标系下的空间分辨率不均匀的雷达数据插值到统一的笛卡尔坐标下形成空间分辨率均匀的互不相关的网格点数据；

第三，利用OPENMP程序将得到的网格点数据在各自线程的控制下进入对应的CPU核中运行即可得到相应的气象数据；

最后，将各CPU核计算所得数据进行合成即可得到所需的气象数据。

对球面坐标系数据进行格点化时，首先利用笛卡尔坐标系下的网格点的高度、格点位置计算出球坐标系中对应仰角、方位和斜距，然后根据计算出的仰角、方位、斜距在雷达球坐标系中的位置，利用内插方法给出该网格点的赋值，得到该网格点上的分析值；然后将球坐标系下的雷达数据插值到笛卡尔坐标系下的细网格上时，用径向和方位上的最近邻居发和垂直线性内插法(nearest neighbor on range-azimuth planes combined with a linearinterpolation in vertieal direction，简称NVI)方法即可得到空间连续的格点数据。

所述的NVI方法是：设定(r，a，e)是某一网格点在雷达球坐标系中的位置，r为斜距，a为方位角，e为仰角。e位于其上下相邻仰角e1和e2之间。(r，a，e₁)和(r，a，e₂)分别是经过该网格点的垂线(仰角低于20°时，垂直方向可用仰角方向近似)与其上下仰角波束轴线的交点，那么该网格点的分析值f^a(r，a，e)可以用这两点的分析f^a(r，a，e₁)和f^a(r，a，e₂)进行垂直线性内插得到，即：

f^a(r，a，e₁)＝[w_e1f^a(r，a，e₁)+w_e2f^a(r，a，e₂)]·(w_e1+w_e2)^-1

其中w_e1和w_e2为内插权重：

w_e1＝(e₂-e)/(e₂-e₁)

w_e2＝(e-e₁)/(e₂-e₁)

f^a(r，a，e₁)和f^a(r，a，e₂)为最靠近点(r，a，e₁)和(r，a，e₂)的雷达距离库的观测值，它们的获取使用了径向和方位上的最近邻居法。r_i-1，r_i，r_i+1为相邻距离库，a_i-1，a_i，a_i+1为相邻方位角，由雷达扫描时半功率线和半功率库所围成的梯形区是距离库r_i的影响区，在径向、方位方向上落在这个梯形区域的点(r，a)得分析值f^a(r，a)都用距离库r_i的观测值f^o(r_i，a_i)来赋值，即：

f^a(r，a)＝f^o(r_i，a_i)。

本发明的有益效果：

本发明通过共享内存并行OpenMP多核并行技术，实现了对天气雷达回波的单雷达二次格点产品和多雷达拼图格点产品计算和生成，并通过对本算法模块的独立封装，具有极强的可移植性和扩展性。

进行雷达数据格点化计算的时候，在初始化后每一个仰角上的计算量比较小，OpenMP自身的调度会消耗一定的时间，所以在对每个仰角进行计算的时候，由于计算量小，使用OpenMP的优势并不明显，有时候甚至由于调度和CPU被其他进程占用的原因，导致并行化后的效率不升反降，例如第6个仰角和第8个仰角，如表1所示。而初始化的时候，数据的计算量约是九个仰角计算量的和，并行化调度的时间占用整个计算过程的时间的比重很小，此时，OpenMP在双核CPU上的计算优势则凸现了出来。

表1OpenMP并行化算法与常规算法的对比

可以看出，在双核的CPU上，通过多核计算并行化后，格点化产品计算时间明显减少，计算速度明显提高。尤其在数据初始化的时候加速比达到了1.71465。

未使用多核计算之前，在双核CPU上，进行产品计算时，CPU的最大使用率只有50％，两个计算核中的一个计算核处于休眠状态，使用本发明的多核并行技术之后，CPU的最大使用率达到97％，充分使用了计算机的运算能力。因此，本发明的OpenMP多核并行技术在天气雷达算法应用方面具有推广的价值。

附图说明

图1是现代计算机共享内存多处理器体系的结构示意图。

图2是本发明的雷达格点生成时的FORK-JOIN模型。

图3是本发明的VCP21扫描模式扫描示意图。

图4是本发明的径向数据的格点化三维格点化示意图。

图5是本发明的多核并行运算示意图。

具体实施方式

下面结合附图和实施例对本发明作进一步的说明。

如图1-5所示。

一种天气雷达数据的多核并行计算方法，它包括以下步骤：

首先，将探测到的天气雷达观测数据用球面坐标系进行存储；

其次，将存储的球面坐标系数据进行格点化处理，也就是把球面坐标系下的空间分辨率不均匀的雷达数据插值到统一的笛卡尔坐标下形成空间分辨率均匀的互不相关的网格点数据；

第三，利用OPENMP程序将得到的网格点数据在各自线程的控制下进入对应的CPU核中运行即可得到相应的气象数据；

最后，将各CPU核计算所得数据进行合成即可得到所需的气象数据。

下面以偏振雷达为例详述如下：

偏振雷达观测资料(每个数据文件)中包含下列几种基本要素：

1)反射率因子：ZH；

2)差分反射率因子：ZDR；

3)差分相位：φDP；

4)径向速度：Vel；

5)速度谱宽：Spw；

6)相关系数：ρHV；

7)差分传播常数：KDP(根据φDP得到)；

双线偏振雷达具有极大的空间分辨率，雷达在进行扫描时，距离分辨率为150M，最大探测距离为150KM。因此，每根径向上每种类型的数据的个数为1000个，一个PPI扫描可以产生365个径向，一个体扫中有9层的PPI数据，因此一个体扫结束，所产生的数据的个数约为：7*1000*365*9＝22995000，各种天气雷达的单雷达二次格点产品和多雷达的拼图产品对应的算法需要基于这些数据进行各种各样的数学计算，由于多核CPU的发展，因此发明了一种多核并行的雷达格点数据产品的计算方法，以充分利用CPU的计算能力来完成天气雷达产品的计算。

雷达做立体观测时，雷达波速是在各个仰角沿方位进行径向扫描，因此，原始雷达观测资料常常用球坐标系(仰角、方位角和斜距)进行存储，而在综合分析雷达资料和其他观测资料(例如卫星)或大范围的拼接多个雷达的资料时，使用球坐标系存储十分不方便。因此必须把球坐标的全分辨率的雷达资料内插到笛卡尔坐标系下的三维格点上。

天气雷达的实际目标是在无人值守的情况下，可以不间断的连续工作，而体扫模式是最佳的扫描模式，雷达以体扫模式运行时，有4种体扫覆盖模式：VCP11、VCP21、VCP31和VCP32。

以标准大气折射下VCP21扫描模式的波束传播路径来说明雷达数据的扫描方式，雷达扫描的示意图如图3所示，图3中虚线为波束轴线，可见雷达观测资料的空间分辨率很不均匀，在方位上的分辨率随着斜距线性增加，同一高度相邻仰角的数据的水平间隔和高度和仰角有很大关系。

雷达资料的格点化时，关键点就是要把球面坐标系下的空间分辨率不均匀的雷达资料插值到统一的笛卡尔坐标下形成空间分辨率均匀的网格点资料，径向数据的的网格化如图4所示。在数据内插中应尽可能保留原始体扫资料中原有的结构特征。

我们利用笛卡尔坐标系下的网格点的高度，格点位置计算出球坐标系中的仰角、方位和斜距，然后根据计算出的仰角、方位、斜距在雷达球坐标系中的位置，利用内插方法给出该网格点的赋值，得到该网格点上的分析值。

把球坐标系下的雷达资料插值到笛卡尔坐标系下的细网格上时，用NVI(径向和方位上的最近邻居发和垂直线性内插法)方法技能得到空间比较连续的格点数据，同时也能最好地保留体扫资料中的原有反射率结构特征.

NVI算法的思想如下：

设定(r，a，e)是某一网格点在雷达球坐标系中的位置，r为斜距，a为方位角，e为仰角。e位于其上下相邻仰角e1和e2之间。(r，a，e₁)和(r，a，e₂)分别是经过该网格点的垂线(仰角低于20°时，垂直方向可用仰角方向近似)与其上下仰角波束轴线的交点，那么该网格点的分析值f^a(r，a，e)可以用这两点的分析值f^a(r，a，e₁)和f^a(r，a，e₂)进行垂直线性内插得到，即：

f^a(r，a，e₁)＝[w_e1f^a(r，a，e₁)+w_e2f^a(r，a，e₂)]·(w_e1+w_e2)^-1

其中w_e1和w_e2为内插权重：

w_e1＝(e₂-e)/(e₂-e₁)

w_e2＝(e-e₁)/(e₂-e₁)

f^a(r，a，e₁)和f^a(r，a，e₂)为最靠近点(r，a，e₁)和(r，a，e₂)的雷达距离库的观测值，它们的获取使用了径向和方位上的最近邻居法。r_i-1，r_i，r_i+1为相邻距离库，a_i-1，a_i，a_i+1为相邻方位角，由雷达扫描时半功率线和半功率库所围成的梯形区是距离库r_i的影响区，在径向、方位方向上落在这个梯形区域的点(r，a)得分析值f^a(r，a)都用距离库r_i的观测值f^o(r_i，a_i)来赋值，即：

f^a(r a)＝f^o(r_i，a_i)

计算过程中，任意格点数据之间的计算互不相关，因此符合并行化的要求，算法可以采用并行的方式实现。雷达格点产品生成时的Fork-Join模型如图2所示。具体的计算过程如图3所示。与图3相配的多行并行运算原理如图1所示，

为了更好的说明，本文通过与常规模式下的回波强度三维格点化计算时间对比，其中的格点数据为某一高度的，计算时间是从格点初始化初始化到9层PPI数据全部计算完毕。程序在多核CPU上执行时，运行的方式如图3所示，线程的个数根据CPU核的个数动态生成，每个线程占用各自的CPU核进行计算。

由于格点数会影响格点化数据产品的计算生成速度，为了充分利用雷达数据资料，我们采用能够达到的最大分辨率，因此设定的格点数为2000×2000，空间范围为300KM×300KM，空间分辨率为150M×150M。

本发明未涉及部分均与现有技术相同或可采用现有技术加以实现。

Claims (3)

1.一种天气雷达数据的多核并行计算方法，其特征是它包括以下步骤：

首先，将探测到的天气雷达观测数据用球面坐标系进行存储；

其次，将存储的球面坐标系数据进行格点化处理，也就是把球面坐标系下的空间分辨率不均匀的雷达数据插值到统一的笛卡尔坐标下形成空间分辨率均匀的互不相关的网格点数据；

第三，利用OPENMP程序将得到的网格点数据在各自线程的控制下进入对应的CPU核中运行即可得到相应的气象数据；

最后，将各CPU核计算所得数据进行合成即可得到所需的气象数据。

2.根据权利要求1所述的天气雷达数据的多核并行计算方法，其特征是对球面坐标系数据进行格点化时，首先利用笛卡尔坐标系下的网格点的高度、格点位置计算出球坐标系中对应仰角、方位和斜距，然后根据计算出的仰角、方位、斜距在雷达球坐标系中的位置，利用内插方法给出该网格点的赋值，得到该网格点上的分析值；然后将球坐标系下的雷达数据插值到笛卡尔坐标系下的细网格上时，用NVI方法即可得到空间连续的格点数据。

3.根据权利要求2所述的天气雷达数据的多核并行计算方法，其特征是所述的NVI方法是：设定(r，a，e)是某一网格点在雷达球坐标系中的位置，r为斜距，a为方位角，e为仰角。e位于其上下相邻仰角e1和e2之间。(r，a，e₁)和(r，a，e₂)分别是经过该网格点的垂线(仰角低于20°时，垂直方向可用仰角方向近似)与其上下仰角波束轴线的交点，那么该网格点的分析值f^a(r，a，e)可以用这两点的分析f^a(r，a，e₁)和f^a(r，a，e₂)进行垂直线性内插得到，即：

f^a(r，a，e₁)＝[w_e1f^a(r，a，e₁)+w_e2f^a(r，a，e₂)]·(w_e1+w_e2)^-1

其中w_e1和w_e2为内插权重：

w_e1＝(e₂-e)/(e₂-e₁)

w_e2＝(e-e₁)/(e₂-e₁)

f^a(r，a，e₁)和f^a(r，a，e₂)为最靠近点(r，a，e₁)和(r，a，e₂)的雷达距离库的观测值，它们的获取使用了径向和方位上的最近邻居法。r_i-1，r_i，r_i+1为相邻距离库，a_i-1，a_i，a_i+1为相邻方位角，由雷达扫描时半功率线和半功率库所围成的梯形区是距离库r_i的影响区，在径向、方位方向上落在这个梯形区域的点(r，a)得分析值f^a(r，a)都用距离库r_i的观测值f^o(r_i，a_i)来赋值，即：

f^a(r，a)＝f^o(r_i，a_i)。

Images