謝千里，林 光，朱偉偉

（中國(guó)民用航空中南地區(qū)空中交通管理局，廣東 廣州510403）

0 引言

天氣雷達(dá)對(duì)于氣象預(yù)報(bào)的作用極其重要，不可替代。隨著我國(guó)天氣雷達(dá)組網(wǎng)建設(shè)，全國(guó)大部分地區(qū)基本完成覆蓋，目前全國(guó)氣象局的新一代天氣雷達(dá)數(shù)量已達(dá)二百多部，民航空管和機(jī)場(chǎng)的天氣雷達(dá)也有一百多部，對(duì)如此多數(shù)量的雷達(dá)數(shù)據(jù)進(jìn)行處理，對(duì)計(jì)算機(jī)的運(yùn)算性能要求越來(lái)越高。CPU擅長(zhǎng)邏輯控制，串行的運(yùn)算和具有復(fù)雜計(jì)算步驟的處理。而GPU則擅長(zhǎng)大吞吐量、大規(guī)模、簡(jiǎn)單運(yùn)算的并發(fā)計(jì)算。天氣雷達(dá)基數(shù)據(jù)解析處理的特點(diǎn)是數(shù)據(jù)量大，計(jì)算密集，重復(fù)量大，但是邏輯計(jì)算簡(jiǎn)單，所以非常適合用GPU進(jìn)行并行計(jì)算。本文提出一種使用GPU來(lái)代替?zhèn)鹘y(tǒng)的CPU進(jìn)行天氣雷達(dá)基數(shù)據(jù)的解析運(yùn)算的方法，大大提高了數(shù)據(jù)的處理速度。

1 天氣雷達(dá)體掃數(shù)據(jù)的結(jié)構(gòu)

如1所示，天氣雷達(dá)體掃數(shù)據(jù)一般由5-14個(gè)仰角層組成，每層一般由360條由中心向四周輻射的徑向數(shù)據(jù)條組成，每條徑向數(shù)據(jù)由從中心到邊沿順序排序的基本單元數(shù)據(jù)組成。體掃基數(shù)據(jù)所覆蓋的空間是一個(gè)很扁的倒圓錐體，但是數(shù)據(jù)是離散的，在這個(gè)錐體中存在大量數(shù)據(jù)空隙。并且數(shù)據(jù)以極坐標(biāo)方式排列，對(duì)這些數(shù)據(jù)進(jìn)行統(tǒng)計(jì)和圖形化顯示，必須變換為直角坐標(biāo)系，并對(duì)空隙進(jìn)行插值填充，以得到一個(gè)柵格狀的類(lèi)似魔方的立體數(shù)據(jù)。

圖1 天氣雷達(dá)體掃數(shù)據(jù)結(jié)構(gòu)

顯然，數(shù)據(jù)中所有點(diǎn)的計(jì)算都可以獨(dú)立進(jìn)行，并不依賴(lài)于其它點(diǎn)處理的結(jié)果。故其運(yùn)算適合進(jìn)行并行分解。

2 英偉達(dá)G P U的硬件結(jié)構(gòu)

如圖2所示，一個(gè)GPU由若干個(gè)SM（流多處理器）構(gòu)成，每個(gè)SM由若干個(gè)SP（流處理器，也稱(chēng)核心）組成，每個(gè)SM中的SP共用一個(gè)寄存器文件（高速寄存器），以及一個(gè)共享內(nèi)存，SP由線程束調(diào)度器進(jìn)行線程調(diào)度，全局內(nèi)存（即顯存）為所有SP所共用。

圖2 英偉達(dá)G P U的硬件架構(gòu)

3 英偉達(dá)C U D A程序結(jié)構(gòu)

CUDA的操作包含3個(gè)基本步驟[1]：

（1）CPU在計(jì)算機(jī)內(nèi)存上準(zhǔn)備好數(shù)據(jù)，在GPU上分配顯存，并把數(shù)據(jù)從內(nèi)存?zhèn)魉偷紾PU顯存。

（2）CPU啟動(dòng)GPU核函數(shù)，GPU創(chuàng)建并發(fā)線程，運(yùn)行核函數(shù)，處理顯存數(shù)據(jù)。

（3）CPU把數(shù)據(jù)從GPU顯存取回到內(nèi)存中，并釋放GPU顯存。

圖3 英偉達(dá)C U D A線程結(jié)構(gòu)

CUDA程序分為host和device兩部分，結(jié)構(gòu)如圖3所示，host運(yùn)行于CPU端，device運(yùn)行于GPU端。CPU通過(guò)核函數(shù)（kernel）啟動(dòng)GPU多線程運(yùn)算。每次核函數(shù)的調(diào)用使用一個(gè)線程網(wǎng)格（grid），一個(gè)線程網(wǎng)格由若干個(gè)線程塊（block）組成，每個(gè)線程塊由若干個(gè)線程（thread）組成。

4 計(jì)算任務(wù)

如圖4所示，計(jì)算任務(wù)是要把有很多空隙的倒圓錐體掃極坐標(biāo)系數(shù)據(jù)通過(guò)插值和變換，得到一個(gè)密實(shí)的柵格化的立方體。

圖4 天氣雷達(dá)體掃數(shù)據(jù)柵格化

計(jì)算每一個(gè)點(diǎn)的取值，該點(diǎn)的取值需要從兩個(gè)方向進(jìn)行計(jì)算，先計(jì)算出該點(diǎn)在相鄰的仰角掃描層表面上的等半徑投影位置，如圖5，然后將上一步得到的每個(gè)投影點(diǎn)位置在該仰角層中再向相鄰兩條徑向線進(jìn)行等半徑投影，得到兩個(gè)徑向條上的點(diǎn)，如圖6[2]。

圖5 垂直方向柵格插值

圖6 水平方向徑向定位

如果計(jì)算點(diǎn)處于兩個(gè)仰角層之間，那么通過(guò)上述投影取值，將獲得4個(gè)最終的投影點(diǎn)位置和值；如果計(jì)算點(diǎn)處于最低仰角之下，或處于最高仰角之上，那么通過(guò)上述投影取值，將獲得2個(gè)最終的投影點(diǎn)位置和值。最后通過(guò)距離加權(quán)計(jì)算出插值的大小。

確定插值后，還要確定一個(gè)標(biāo)準(zhǔn)的位置，即進(jìn)行海拔高度和水平距離的計(jì)算。如圖7所示，O點(diǎn)為雷達(dá)天線所在位置，h0為雷達(dá)天線所處的海拔高度，A點(diǎn)為探測(cè)范圍空中某點(diǎn)，B點(diǎn)為A點(diǎn)與地心的連線和海平面的交點(diǎn)。r為A點(diǎn)距離雷達(dá)天線的直線距離，R為地球半徑，線段AB即為A點(diǎn)的海拔高度[2]。

圖7 柵格海拔高度計(jì)算

A點(diǎn)所處位置的海拔高度AB由兩個(gè)部分組成：

式中，h1為A點(diǎn)距雷達(dá)天線所處的海平面延伸平面的垂直高度；h2為因地球曲率而增加的高度；A點(diǎn)距雷達(dá)站的水平距離為弧BC的長(zhǎng)度B（C。

由幾何關(guān)系推導(dǎo)出：

由（1）（2）（3）（4）式可以計(jì)算出某點(diǎn)的海拔高度和離雷達(dá)站的水平距離，通過(guò)該算法進(jìn)行三維柵格化[2]。

5 計(jì)算任務(wù)的并行分解

下面以一個(gè)范圍為150 km，庫(kù)長(zhǎng)為300 m的雷達(dá)，使用英偉達(dá)Quadro M2000M顯卡處理為例，來(lái)設(shè)計(jì)一個(gè)并行算法。

柵格化后的立方體每一小格的尺寸為300 m邊長(zhǎng)的正方體，高度取15 000 m，則該立方體的長(zhǎng)寬高為：1 000× 1 000× 50，共 50，000，000 （約 50 M）個(gè)數(shù)據(jù)值。

每個(gè)柵格點(diǎn)的取值計(jì)算由一個(gè)線程運(yùn)行。每個(gè)線程計(jì)算出某一個(gè)位置的柵格點(diǎn)的取值。算法見(jiàn)第4部分計(jì)算任務(wù)的描述。

GPU的基本調(diào)度單元是線程束，一個(gè)線程束需要占用一個(gè)SM來(lái)運(yùn)行，多個(gè)線程束需要輪流進(jìn)入SM，對(duì)于M2000M，SM的數(shù)量是5，線程束大小是32，即一個(gè)線程束由32個(gè)線程組成，故每個(gè)線程塊中線程的數(shù)量應(yīng)設(shè)計(jì)為32的整數(shù)倍，才能發(fā)揮最大的并發(fā)性。M2000M的硬并發(fā)線程數(shù)量實(shí)際上是5×32=160個(gè)。每個(gè)線程塊有獨(dú)立的寄存器和共享內(nèi)存，供一個(gè)塊中的線程共同使用，故每個(gè)線程塊中的線程數(shù)量不宜過(guò)多。例如對(duì)于M2000M，每個(gè)線程塊最大線程數(shù)為1024，那么，對(duì)于本計(jì)算任務(wù)，每個(gè)線程塊中的線程數(shù)量設(shè)為128個(gè)較為合理。線程的維數(shù)和線程塊的維數(shù)均為邏輯劃分，對(duì)于實(shí)際性能并無(wú)影響，故為了簡(jiǎn)化數(shù)據(jù)結(jié)構(gòu)以便于向GPU傳送和取回?cái)?shù)據(jù)，把數(shù)據(jù)1維化，所以線程塊和線程的組織均只需1個(gè)維度即可。

那么核函數(shù)kernelGet3dVolume<<>>（…）中，塊中線程數(shù)量threadNum的值為128，則線程塊數(shù)量blockNum=50，000，00/128=390625。

6 核函數(shù)的設(shè)計(jì)

核函數(shù)原型如下：

__global__void kernelGet3dVolume（int layer－Count， int radialCount， int stdBinCount， int std－HeightCount， int stdBinWidth， int altitude， float stdSweepAngle，float*pmtElevation，int*z，unsigned char*volume，int volumeTotalCount）

參數(shù)說(shuō)明：

layerCount：體掃數(shù)據(jù)層數(shù)

radialCount：每層徑向條數(shù)

stdBinCount：每條徑向庫(kù)數(shù)

stdHeightCount：柵格化后的z軸柵格數(shù)

stdBinWidth：庫(kù)長(zhǎng)

altitude：雷達(dá)站海拔高度

stdSweepAngle：水平掃描角度間隔

pmtElevation：各層仰角

z：按順序排列的庫(kù)數(shù)據(jù)，從低仰角到高仰角，正北徑向開(kāi)始順時(shí)針旋轉(zhuǎn)，在徑向中從最近點(diǎn)開(kāi)始到最遠(yuǎn)點(diǎn)。

volume：處理完成后生成的三維柵格數(shù)據(jù)

volumeTotalCount：volume的字節(jié)數(shù)。

該核函數(shù)的作用是計(jì)算三維柵格中某一個(gè)點(diǎn)應(yīng)取的值。

回波的值在z中，z是一維數(shù)組，可以通過(guò)體掃數(shù)據(jù)層數(shù)、每層徑向條數(shù)等其它參數(shù)來(lái)配合數(shù)組下標(biāo)計(jì)算，定位出在倒錐面中的空間位置，以進(jìn)行插值計(jì)算。

CPU整理好體掃數(shù)據(jù)后裝載至一塊連續(xù)內(nèi)存中，然后將該內(nèi)存塊復(fù)制至GPU內(nèi)存即顯存中的數(shù)組z，啟動(dòng)核函數(shù)展開(kāi)所有線程，產(chǎn)生50，000，000個(gè)邏輯并發(fā)線程，由GPU按blockNum和threadNum兩個(gè)參數(shù)來(lái)進(jìn)行自動(dòng)調(diào)度。

7 G P U與C P U的性能對(duì)比

下面對(duì)714CDP雷達(dá)基數(shù)據(jù)201607061059240.05V進(jìn)行處理，文件信息如圖5所示。

圖5 714CD P型天氣雷達(dá)基數(shù)據(jù)

該體掃總層數(shù)為14，每層徑向數(shù)為720，每條徑向的庫(kù)數(shù)為1 000。柵格化后生成組合反射率產(chǎn)品如圖8所示。

圖8 基數(shù)據(jù)201607061059240.05V生成的組合反射率圖

CPU和GPU均使用C++代碼，CPU代碼已進(jìn)行多線程優(yōu)化，已充分利用CPU的多核，編譯器均為VS2017，對(duì)某個(gè)天氣雷達(dá)進(jìn)行1 000×1 000×50個(gè)點(diǎn)的柵格化計(jì)算，使用相同的處理算法和天氣雷達(dá)基數(shù)據(jù)文件，分別在筆記本電腦（聯(lián)想Thinkpad P50），以及臺(tái)式機(jī)（聯(lián)想Thinkstation P510）上多次運(yùn)行，得到耗時(shí)平均值如表1所示。

表1 G P U和C P U處理天氣雷達(dá)基數(shù)據(jù)性能實(shí)測(cè)

在筆記本電腦上，E3-1505是移動(dòng)端的高端CPU，M2000M是移動(dòng)端的專(zhuān)業(yè)繪圖顯卡。在臺(tái)式機(jī)上，E5-1620是工作站級(jí)別的中高端性能 CPU，P2000是中端繪圖顯卡，故它們之間的對(duì)比是匹配的，具有參考價(jià)值的。

從上表可以看出，在筆記本和臺(tái)式機(jī)上，GPU的運(yùn)算速度均比CPU高約8.5倍。假如有200部天氣雷達(dá)體掃基數(shù)據(jù)匯總到一臺(tái)主機(jī)上處理，那么，CPU所需時(shí)間為1.038×200=208 s，即3分28秒，而使用GPU，則僅需0.122×200=24 s，比CPU節(jié)約3分4秒。

8 結(jié)束語(yǔ)

GPU用于高性能計(jì)算已有多年，其海量硬并發(fā)的特點(diǎn)，很適合用于天氣雷達(dá)基數(shù)據(jù)的處理，比起用CPU計(jì)算，速度得到顯著的提升，如果要處理的天氣雷達(dá)數(shù)據(jù)很多，例如200部，那么，使用GPU后，所累積使用時(shí)間將大大縮減，即氣象產(chǎn)品的處理延時(shí)大大縮減，使氣象預(yù)報(bào)工作能更及時(shí)地開(kāi)展。