天氣雷達體掃數據的G P U并行處理

謝千里，林 光，朱偉偉

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

0 引言

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

1 天氣雷達體掃數據的結構

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

圖1 天氣雷達體掃數據結構

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

2 英偉達G P U的硬件結構

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

圖2 英偉達G P U的硬件架構

3 英偉達C U D A程序結構

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

（1）CPU在計算機內存上準備好數據，在GPU上分配顯存，并把數據從內存傳送到GPU顯存。

（2）CPU啟動GPU核函數，GPU創建并發線程，運行核函數，處理顯存數據。

（3）CPU把數據從GPU顯存取回到內存中，并釋放GPU顯存。

圖3 英偉達C U D A線程結構

CUDA程序分為host和device兩部分，結構如圖3所示，host運行于CPU端，device運行于GPU端。CPU通過核函數（kernel）啟動GPU多線程運算。每次核函數的調用使用一個線程網格（grid），一個線程網格由若干個線程塊（block）組成，每個線程塊由若干個線程（thread）組成。

4 計算任務

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

圖4 天氣雷達體掃數據柵格化

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

圖5 垂直方向柵格插值

圖6 水平方向徑向定位

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

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

圖7 柵格海拔高度計算

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

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

由幾何關系推導出：

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

5 計算任務的并行分解

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

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

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

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

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

6 核函數的設計

核函數原型如下：

__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）

參數說明：

layerCount：體掃數據層數

radialCount：每層徑向條數

stdBinCount：每條徑向庫數

stdHeightCount：柵格化后的z軸柵格數

stdBinWidth：庫長

altitude：雷達站海拔高度

stdSweepAngle：水平掃描角度間隔

pmtElevation：各層仰角

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

volume：處理完成后生成的三維柵格數據

volumeTotalCount：volume的字節數。

該核函數的作用是計算三維柵格中某一個點應取的值。

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

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

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

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

圖5 714CD P型天氣雷達基數據

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

圖8 基數據201607061059240.05V生成的組合反射率圖

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

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

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

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

8 結束語

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