傅 游,徐方正,梁建國
(山東科技大學 計算機科學與工程學院,山東 青島 266590)
平均海平面(mean sea surface,MSS)是指相對于參考橢球在一定時間內的平均動態海面高,由平均海面地形和大地水準面兩部分相加得到[1]。MSS在研究地殼形變、大洋環流、海洋重力計算[2]、大地水準面起伏確定和地殼形變[3]等問題中得到廣泛應用,對地球科學和環境科學研究具有重要意義。
建立MSS模型面臨的挑戰是在有限的時間跨度內實現對時間海面變化的最精確濾波,同時獲得最高的空間分辨率。通常結合來自精確重復任務(exact repeat mission,ERM)的數據與ERS-1、GEOSAT等早期的大地測量任務(geodetic mission,GM)測高數據實現[1]。校正GM數據中的海面時變信號時,需要讀取大量的衛星軌跡數據文件[3],比如在構建日本海周邊區域的平均海平面模型中對GM數據進行海面時變校正時,共使用了5顆衛星的GM衛星測高數據及其對應的ERM測高數據[4]。以Crystat-2衛星(2011.01.28—2019.12.12)為例,其GM衛星測高數據包含112個周期,總共92 669個軌跡文件。生成結果的精度和分辨率要求越高,I/O的數據越多,計算耗時和I/O耗時越長,必須利用并行計算技術進行加速。并行計算已在衛星數據處理領域發揮了重要作用[5-7]。在I/O性能優化方面,Schenck等[8]提出將突發式數據緩存使用快速存儲介質作為緩沖區,將進程之間由I/O引起的負載不平衡降到最低限度,同時加快I/O的整體速度;Thakur等[9]提出的ROMIO使用雙階段I/O調度算法優化后的集中式I/O,盡可能對同一文件但不同數據段的多次訪問進行合并以減少訪問次數;Behazad等[10]使用遺傳算法考慮影響I/O性能的各個參數,再進行全局空間搜索,從而尋找最優參數解;Chen等[11]使用遺傳算法對并行I/O性能進行自動調優;Guedes等[12]針對運行在基于容器的服務器虛擬化集群上的I/O密集型應用進行研究,在虛擬環境下提供緩存服務,在大規模集群的存儲文件系統中(如Lustre、通用并行文件系統和Panasas文件系統),將單個文件分為多個子文件存儲在多個數據服務器上,通過服務器的并發來提高I/O效率。這些方法在一定程度上確實能減少I/O耗時,但不同程序具有不同I/O特征,要取得更好的I/O優化效果,必須針對具體程序進行分析,制訂具體優化策略。
在全球平均海面模型的建立過程中,在Yuan等[13]研究基礎上,開發了基于時空客觀分析法的GM衛星測高數據海面時變校正的串行程序,完成了向高性能集群系統的移植。該串行程序讀取的多源衛星總數據量約2 TB,輸出數據約500 GB,總讀取軌跡文件數約10 000萬個。在CPU為Intel i7-10875H、內存16 GB的個人電腦上需運行約3個月,嚴重影響研究進度。而在完成高性能集群系統移植后,程序計算時間大大減少,但I/O作為系統性能瓶頸的情況并未改善,影響了系統可擴展性。
為了縮短I/O耗時,實現系統可擴展性,本研究從兩方面對GM衛星測高數據海面時變校正程序I/O特征進行分析;為了提高I/O效率,提出按周期分配方案,并針對該方案可擴展性不佳、易導致負載不均衡的問題,提出一種合并再分配方案;使用消息傳遞接口(massage passing interface,MPI)文件視口函數對合并再分配算法進行優化,進一步提高I/O效率。
以進行海面時變校正的GM衛星測高數據在相同時間跨度的ERM測高數據為參考基準,先將ERM測高數據的海平面異常(sea level anomalies,SLA)與GM衛星測高數據進行時空匹配,再進行海面時變校正,即可得到校正后的GM衛星測高數據。采用時空客觀分析法進行GM衛星測高數據海面時變校正的關鍵是如何選取與待校正GM衛星測高數據在時間和空間相匹配的ERM測高數據并計算SLA。GM衛星測高數據海面時變校正串行算法流程如圖1所示,圖中左側為ERM的SLA數據的篩選和計算過程,右側為GM衛星測高數據的篩選和計算過程。

圖1 串行算法流程圖
在對原串行Fortran程序使用插樁法進行熱點檢測后,發現該串行程序讀寫文件耗時占比達80.44%,而數據計算部分僅占10.37%。若能優化程序的I/O部分,則該程序的整體耗時會大幅減少。
時空客觀分析法的特征是在考慮時空尺度的前提下,將沿測高軌跡的SLA數據格網化為規則的格網SLA,再對同時空范圍內的GM衛星測高數據進行校正。GM和ERM衛星測高數據文件結構的特點是包含多個周期文件夾且每個周期文件夾包含多個軌跡文件,文件目錄結構如圖2所示。

圖2 衛星測高數據文件目錄結構圖
軌跡文件中以文本形式記錄衛星測高數據,GM衛星測高數據和ERM衛星測高數據文件內容結構分別如圖3、圖4所示。圖3和圖4中,每一行代表一個觀測點的信息,每一列代表一個屬性,其他信息篩除。GM衛星測高數據的觀測點信息屬性包括觀測時刻、經度、緯度、動態海面高、大地水準面和平均動態海面高,其中大地水準面和平均動態海面高在本研究的海面時變校正中未被采用,但在后續海面高建模中采用,因此未去除。ERM衛星測高數據包含時刻、經度、緯度和SLA。

圖3 GM衛星測高數據文件內容結構圖

圖4 ERM衛星測高數據文件內容結構圖
GM衛星測高數據海面時變校正程序I/O占比大的主要原因包括:①測高數據分布在數量繁多的軌跡文件中,如引言中提到的Crystat-2衛星包含92 669個軌跡文件,而對應的同時期ERM軌跡文件83 058個,共計175 727個軌跡文件;②讀入GM衛星測高數據和查找時空對應的SLA數據過程中,需要多次讀入不同文件,頻繁切換文件句柄,切換頻率約5 000次/s。每次完成一個觀測點的計算均需要寫入文件,每次寫入的數據量僅48字節,輸出文件時需寫入約10 000次/s??梢奊M測高數據時空客觀分析法程序具有I/O密集型程序的特征,頻繁的I/O導致程序運行速度受限于I/O帶寬,無法充分發揮大規模集群中多核計算機性能。
為了將該程序移植到高性能集群系統,實現多進程并行I/O,本研究提出一種按周期分配數據的并行方案。在該方案中,以一個周期的GM衛星測高數據為最小任務分配粒度,即每個進程處理整數個周期的GM衛星測高數據,盡可能把所有數據按周期數平均分配給每個進程進行計算,兩兩進程之間分配到的周期數差最大為1。設n為總進程數,T為總任務數,進程i處理的任務數為ti,
(1)
式中:/為取商運算,%為取余運算。
在多進程運行時發現,各進程運行時間會有較大差異。以10個不完整周期(包含約100個軌跡文件)的GM衛星測高數據為測試集,該方案并行程序分別開啟1到10個進程的運行時間與加速比如圖5所示,9進程時各進程運行時間如圖6所示。

圖5 并行程序運行時間與加速比圖

圖6 并行程序9進程時各進程運行時間
由圖5可以看出,當進程數由5升至9時,運行時間基本保持不變;只有當進程數為1、2、5和10時,加速比與進程數接近。由圖6可以看出,0#進程的運行時間約為其他進程的2倍,產生了明顯的負載不均衡。由式(1)可知,當T不能被n整除時,分配給各個進程的任務數不均等,任務差為一個周期的GM衛星測高數據校正任務。任務差越大,負載不均衡的情況越明顯。由于GM衛星測高數據各周期的數據量不等,即使每個進程都分配相同數量的周期,實際的任務數據量差別也很大。除此之外,按周期分配方案的可擴展性不佳,能開啟的進程數不能大于周期數。
由1.2節中的數據文件結構分析可知,若將行數據均勻地劃分給每個進程,需將所有文件的行數進行匯總,統計每個進程的行數。由于數據分布在多個軌跡文件中,在進行SLA數據匹配過程中依然存在頻繁切換文件的問題,出現同一個軌跡文件數據被分配給不同進程的情況,降低任務分配的效率。為避免此類問題的發生,本研究提出合并再分配的并行方案。在實際的任務計算中,最小計算任務單位是一個GM衛星測高的觀測點數據,即一行數據,若能以行為進程分配最小任務單位,可大大提高負載均衡程度。
合并再分配的并行方案將所有衛星測高數據和軌跡文件匯總為如圖7所示的natt個屬性文件。圖7中,call_input表示的匯總文件總行數,即各進程需要讀入的行數之和,ci_input表示進程i分到的衛星測高數據行數。將衛星測高數據匯總整合為數個文件,文件個數為需要測高數據中包含的屬性數(如時間、經度、緯度等),并按時間順序排序。ERM數據的匯總屬性文件除了圖4所示的4個屬性外,每行數據還需要加上該觀測點的周期號和軌跡號,以便于時空匹配。其中,

圖7 合并再分配方案圖

(2)
設進程i在匯總文件開始讀數據的起始偏移量
(3)
各進程可以自行計算出ci_input和oi_input,無需通信。
基于合并再分配方案的并行GM衛星測高數據海面時變校正具體實現如算法1所示。

算法1 基于合并再分配方案的并行GM衛星測高數據海面時變校正算法輸入:GM衛星測高數據gm_data,ERM衛星的SLA數據sla_data輸出:校正后的GM測高數據result_data1: Begin2: MPI及其常用參數初始化,獲取進程號i,進程數n3: 多進程讀取sla_data的不同軌跡文件,并將其匯總為6個SLA屬性文件4: 多進程讀取gm_data,并將其匯總為6個GM屬性文件5: 根據i和n確定本進程讀取的gm_data在文件中的起始位置與長度6: 每個進程將本進程gm_data時間差在10天內的SLA數據讀入內存7: 前6號進程讀取sla_data并廣播給所有進程,存入6個SLA屬性數組8: end_line_gm_data←start_line_gm_data + my_gm_data_length9: for start_line_gm_data←my_gm_data_offset to end_line_gm_databy max_size do ∥每次循環取最多max_size 個觀測點的GM數據,max_size為一次最多可讀取觀測點數,自行設置10: block_gm_data←read_gm_data_to_array(start_line_gm_data,max_size,end_line_gm_data)∥將本次循環要校正的gm_data讀入內存11: for grid_gm_data in block_gm_data do ∥按順序從當前的block_gm_data中選取1°×1°格網內的gm_data12: grid_sla_data←choose_sla_data(my_sla_data) ∥ 篩選距離網格中心不超過1 000 km的sla數據13: grid_sla_data←filter_sla_data(grid_sla_data) ∥ 進行濾波和重采樣操作14: grid_result_data←objective_analysis(grid_sla_data,grid_gm_data)∥使用時空客觀分析法計算出校正后的GM衛星測高數據15: block_result_data←add_result_data(block_result_data,grid_result_data)∥將當前格網的結果數據加入到當前周期結果數據16: end for17: write_block_result_data(block_result_data)∥將一個block的結果數據寫入磁盤18: end for19: 0#進程將所有block的結果文件匯總為一個文件并輸出20: End
為了減少合并再分配時對文件匯總的耗時,利用MPI文件視口函數實現多進程I/O并行加速。MPI文件視口會給每一個進程定義一個獨立文件指針,讀寫位置由當前文件指針確定[14],各進程可以同時讀/寫一個文件,互不干擾。
除了讀/寫的文件路徑之外,每個進程在讀/寫前還需要獲取自己在文件中進行讀/寫操作的偏移量和讀/寫的數據量。每個進程讀操作的數據量依據進程數平均分配,并行讀操作的文件起始偏移量可根據其進程號用式(3)計算得到。但是在并行寫操作時,輸出文件中的數據順序需要與進程順序保持一致,而每個進程寫入輸出文件的數據量并不相等,所以進程號大于0的進程進行寫操作時的偏移量需由上一個進程傳入。
各進程偏移量通信流程如圖8所示。圖8中,ci_output為i號進程需要向輸出文件寫入的數據行數,oi_output為i號進程在輸出文件中的起始偏移量,

圖8 寫文件偏移量通信流程
(4)
所有進程的偏移量傳輸過程按進程號從小到大依序進行,在獲得其輸出文件的偏移量后,通過MPI文件視口輸出函數將結果數據并行寫入輸出文件。
傳統的主從并行方案(以下簡稱傳統方案)中,往往只將計算任務并行化,輸入數據時需要由主進程讀文件并向從進程分發數據;輸出數據時需要從進程向主進程發送數據,再由主進程寫文件。并行I/O方法可避免主從進程間進行大量數據通信,理論上將讀/寫時間縮短為傳統方案的1/n。
實驗平臺為17臺多核計算機組成的集群系統,其中16臺為計算節點,1臺為管理節點。計算節點的整體配置相同,均包含2個14核Intel Xeon Gold 6132處理器和96 GB的DDR4內存。因實驗平臺使用多處理器計算(multi-processor computing,MPC)技術,每個核心可開啟一個進程,可開最大進程數為16×28=448。實驗平臺軟件環境如表1所示。

表1 實驗平臺軟件環境
實驗數據集來源于AVISO(archiving,validation and interpretation of satellite oceanographic data)發布的沿測高軌跡的Level-2+(L2P)SLA數據產品(以下簡稱L2P)。L2P包含多個衛星測高數據,實驗使用其中的T/P、Jason-2和Cryosat-2測高衛星數據集。L2P中每個測高衛星均包含若干周期數據,每個周期數據中又包含若干軌跡數據,而每一個軌跡數據均包含了觀測時間、緯度、經度、衛星軌道高度、衛星到星下點距離、SLA、平均海面高、有效數據標志以及各項誤差改正項等。實驗數據集對于GM數據,保留了觀測時刻、經度、緯度、動態海面高、大地水準面和平均動態海面高;對于ERM數據,保留了觀測時刻、經度、緯度和SLA,將該數據集稱為大數據集。小數據量數據集系從L2P中按固定時間間隔抽選的文件得到。兩組數據集的具體信息如表2所示。

表2 衛星測高數據信息表
1) 負載均衡度測試
為了驗證合并再分配方案的負載均衡度的優勢,使用小數據量數據集,收集了在單節點上開啟6~10進程時,按周期分配方案和合并再分配方案的各進程空閑時間之和,結果如圖9所示。

圖9 多進程下不同并行方案進程空閑時間和
由圖9可以看出,按周期分配方案的進程空閑時間之和,在6到9進程時隨著進程數的上升而增加,而在10進程時大幅降低。這是因為10進程時,周期數剛好能被進程數整除。而合并再分配方案的進程空閑時間在6到9進程時遠低于按周期分配方案,只有在10進程時相近,且隨著周期數上升幅度很小,可見合并再分配方案的負載均衡度優于按周期分配方案。
2) 兩種方案和I/O并行加速效果測試
使用原版串行程序利用小數據量數據集進行實驗,耗時為1 027.89 s。
為驗證兩種并行方案以及I/O并行加速的優化效果,使用小數據集分別進行測試。小數據集包含10個周期,使用按周期分配方案時,若進程數超過10會出現閑置進程,故最大進程數不超過10。合并再分配方案中待匹配的ERM數據有6個屬性(時間、經度、緯度、SLA、所在周期、所在軌跡),每個進程讀入一個屬性,故進程數至少為6。測試結果如圖10所示。

圖10 多進程下不同并行方案運行時間
由圖10可以看出,對比串行程序,按周期分配方案在單進程時耗時更長,是因為啟用MPI需要耗時并且多了任務分配的步驟。6到9進程時,按周期分配方案的耗時基本不變,這是因為出現了任務分配不均,只分配到1個周期的進程需要等待2個周期的進程執行結束。
合并再分配方案在6和10進程時比按周期分配方案耗時更長,是因為分配任務的步驟更加復雜,當按周期分配方案中周期數可以被進程數整除或者余數很小時,按周期分配方案負載較為均衡,合并再分配方案的負載均衡優勢沒有體現出來。當進程數為7到9時,按周期分配方案的負載產生了很大的不均衡,耗時大于合并再分配方案。同時,I/O并行加速對于按周期分配方案也有一定的加速效果,在6進程時加速最多,縮短耗時32.47 s,加速13.4%。
3) 強、弱可擴展性測試
為了驗證經I/O并行加速后的合并再分配方案的強、弱可擴展性,分別使用小數據量和大數據量數據集進行測試。
許多并行計算平臺(例如高性能計算集群)能更高效地處理進程數為2的次方數的并行程序,這是因為此類平臺通常使用硬件結構來組織二叉樹或超立方體拓撲中的節點,使用2的次方數作為進程數允許進程自然映射到硬件拓撲中的節點,可以減少通信延遲并提高性能。因合并再分配方案的進程數不得小于6,故最小進程數設置為8。
小數據量實驗組的加速效果如表3所示。表3中,加速比1由串行程序運行時間除以測試算例運行時間得到,加速比2由8進程并行程序運行時間除以測試算例運行時間得到,并行效率由加速比1除以進程數得到。由表3可知,運行時間在64進程時最短,為86.41 s,且加速比1為11.90,加速比2為1.91,但是并行效率最差,為0.19。并行效率在進程數為8時最高,為0.78。

表3 小數據量實驗組加速效果
小數據量實驗組中,在進程數達到16后,運行時間再無明顯縮減,并行效率隨著進程數的增加而減少。8進程到16進程加速比1和加速比2上升明顯,但是進程數繼續增加后,沒有明顯提升。原因是小數據量實驗組的數據量較小,符合Amdahl定律,小數據量實驗不具有良好的強可擴展性,難以體現出并行計算的優勢。
大數據實驗組中數據量大,最小進程數若設置太小會導致運行時間過長,并且可設置最大進程數為448,故設置并行運行的最小進程數為448的1/8,即56。表4為大數據量實驗組的加速效果,加速比由并行程序56進程運行時間除以測試算例運行時間得到。

表4 大數據量實驗組加速效果
由表4可知,運行時間隨著進程數的每次翻倍,加速比的增長速度也接近翻倍,在進程數為448時最短,為9 283.84 s,加速比也在此時達到最高7.21,并且加速比的上升速率沒有明顯降低,說明該算法具有良好的強可擴展性。
分別以8和56為基準進程數,進程數增加為8倍后,大數據量實驗組的加速比提升為小數據量實驗組的3.78倍,這主要是因為數據量擴大后合并再分配數據的耗時占比更小了,可見該并行方法的加速比在隨著數據量和進程數的增加而提升,也具有良好的弱可擴展性。
本研究實現了時空客觀分析法對GM衛星測高數據的海面時變校正的串行程序,分析了I/O密集型程序特性。使用按周期分配的并行方案對其進行并行化,提出合并再分配方案以保證負載均衡,并使用MPI文件視口函數進行I/O并行優化。實驗結果表明:與按周期分配方案相比,合并再分配方案在多進程運行時耗時更少,并且在周期數不能被進程數整除時,負載均衡度更高,多進程可擴展性也更好;I/O并行加速可縮短合并再分配方案的運行時間;I/O并行加速后的合并再分配方案具有良好的強、弱可擴展性。后續工作可以進一步優化程序的協方差矩陣計算步驟,減少計算環節耗時。