天氣雷達體掃數(shù)據(jù)的G P U并行處理

謝千里，林 光，朱偉偉

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 計算任務(wù)

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

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

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

圖5 垂直方向柵格插值

圖6 水平方向徑向定位

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

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

圖7 柵格海拔高度計算

A點所處位置的海拔高度AB由兩個部分組成：

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

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

由（1）（2）（3）（4）式可以計算出某點的海拔高度和離雷達站的水平距離，通過該算法進行三維柵格化[2]。

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

下面以一個范圍為150 km，庫長為300 m的雷達，使用英偉達Quadro M2000M顯卡處理為例，來設(shè)計一個并行算法。

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

每個柵格點的取值計算由一個線程運行。每個線程計算出某一個位置的柵格點的取值。算法見第4部分計算任務(wù)的描述。

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

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

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

核函數(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ù)說明：

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

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

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

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

stdBinWidth：庫長

altitude：雷達站海拔高度

stdSweepAngle：水平掃描角度間隔

pmtElevation：各層仰角

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

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

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

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

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

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

7 G P U與C P U的性能對比

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

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

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

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

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

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

在筆記本電腦上，E3-1505是移動端的高端CPU，M2000M是移動端的專業(yè)繪圖顯卡。在臺式機上，E5-1620是工作站級別的中高端性能 CPU，P2000是中端繪圖顯卡，故它們之間的對比是匹配的，具有參考價值的。

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

8 結(jié)束語

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