基于局地集合变换卡尔曼滤波的全球海洋资料同化系统设计及算法加速
Global Ocean Data Assimilation System Design and Algorithm Acceleration Based on Local Ensemble Transform Kalman Filter
通讯作者: 徐芳华(1979-),女,天津人,副教授,主要从事海洋数值模式、数据同化、海洋后报和预报. E-mail:fxu@tsinghua.edu.cn
收稿日期: 2019-01-21 修回日期: 2019-04-02 网络出版日期: 2019-06-11
基金资助: |
|
Corresponding authors: Xu Fanghua (1979-), female, Tianjin City, Associate professor. Research areas include numerical ocean model, data assimilation, ocean hindcast and forecast. E-mail:fxu@tsinghua.edu.cn
Received: 2019-01-21 Revised: 2019-04-02 Online: 2019-06-11
作者简介 About authors
范峥(1994-),男,安徽蚌埠人,硕士研究生,主要从事地学高性能计算研究.E-mail:fan-z16@mails.tsinghua.edu.cn
通过对局地集合变换卡尔曼滤波(LETKF)算法的计算时间复杂度的完整分析,发现计算集合空间分析场误差协方差的逆矩阵这一过程计算量最大,耗时最长。且在并行计算环境下,该步骤CPU计算量分配不均是影响计算效率的直接原因。为解决这一问题,采用“贪心算法”设计了一套新的负载均衡策略,并使用该策略开发了一个基于LETKF和并行海洋模块2(POP2)的高性能并行海洋资料同化系统。将2004年1~2月日平均的最优插值海表温度资料(OISST)和同时期的Argo温盐剖面资料同化进入POP2。结果表明,同化有效降低了温度和盐度的均方根误差。同时,在不改变计算结果的前提下,相比原始同化系统,新系统计算性能提升1倍。在更高分辨率(0.1°×0.1°)下,该系统的计算性能仍然可以提升1倍,说明新设计的负载均衡方案稳定可靠。该方案具有很强的可扩展性和移植性,在业务预报中有广泛的应用前景。
关键词:
An integrated analysis about computational time complexity of the Local Ensemble Transform Kalman Filter (LETKF) was performed. It is found that the calculation step of inverse matrix of the error covariance in ensemble space is the most computationally intensive and time consuming. In a parallel computing environment, the uneven distribution of CPU calculations in this step directly leads to low computational efficiency. To solve this problem, a new load balancing strategy was designed based on the "greedy algorithm". A high-performance parallel ocean data assimilation system based on the LETKF was developed and tested using this strategy. This system was based on the Parallel Ocean Program 2 (POP2) of the Community Earth System Model (CESM). The optimal interpolated sea surface temperature data (OISST) and Argo temperature profile data from January to February, 2004 were assimilated into the POP2. The results show that data assimilation effectively reduces the root mean square error of temperature and salinity. Using the new strategy, the exact same results are obtained but the computation time is reduced by half. At higher resolution (0.1°×0.1°),the computing performance is still doubled, indicating that this load balancing scheme is stable and reliable. In addition, the new method has high scalability and portability with great potential to be applied in operational forecasting.
Keywords:
本文引用格式
范峥, 李宏, 刘向文, 徐芳华.
Fan Zheng, Li Hong, Liu Xiangwen, Xu Fanghua.
1 引 言
然而,高昂的计算代价始终是制约资料同化技术发展的关键因素,资料同化方法始终在平衡同化效果和计算量之间发展。这主要体现在同化方法和观测资料两方面。在同化方法方面,例如连续资料同化方法中的三维变分(Three Dimensional Variation,3D-Var)与四维变分(Four Dimensional Variation,4D-Var)[3],虽然三维变分计算量小,但由于没有时间变量使得动力模式无法进行约束,加之同化时间窗口内的观测数据都固定为同一时刻,因此同化结果会产生较大偏差[4];而四维变分需要求解伴随模式,计算量大为增加,使其难以大规模应用。而对于顺序类的资料同化方法,通过对方法的改进来降低计算量更是其发展的主线。为使卡尔曼滤波(Kalman Filter,KF)算法[5]能够用于资料同化,众多学者做了大量探索性的工作,以基于蒙特卡洛模拟方法的集合卡尔曼滤波(Ensemble Kalman Filter,EnKF)[6]为代表。Hunt等[7]通过引入局地化与集合变换的策略针对EnKF中矩阵求逆这一计算热点进行优化,提出局地集合变换卡尔曼滤波(Local Ensemble Transform Kalman Filter,LETKF),进一步降低了计算量,并得到了广泛的应用[8,9,10]。在观测资料方面,随着热带大气海洋阵列(Tropical Atmosphere/Ocean array,TAO)[11],地转海洋学实时海洋观测阵列(Array for Real-time Geostrophic Oceanography,Argo)等国际海洋观测计划的进展[12]及各种带有探测海洋环境要素的卫星系统的建立和完善,使观测资料的数量和类型海量增加,时空分布的差异性不断扩大。研究者们一方面想尽可能同化更多的观测资料以改善同化效果,另一方面想在节约计算资源和降低计算成本的同时缩短计算时间,以便提高科学研究特别是业务预测和预报的效率。
以往为解决资料同化计算量过大的问题主要是对方法本身进行改进,虽然可以提高计算性能,但往往会牺牲相当部分的同化效果。与此同时,很少有人仔细分析资料同化系统在并行计算环境下的计算流程和瓶颈,找到影响其计算效率的关键因素并加以解决。基于此,结合LETKF算法非常适用于并行计算的特点[13]、同化效果较好且有较广阔应用前景的事实[14,15],本文选定LETKF算法进行性能优化。我们详细分析了该算法各步骤的计算时间复杂度,发现该算法的瓶颈在于计算集合空间分析场误差协方差的逆矩阵这一步,且CPU计算量分配不均是影响该步骤性能的直接原因。本文首次利用“贪心算法”设计一套负载均衡策略以优化CPU的计算效率,开发了一个基于LETKF的高性能并行海洋资料同化系统。该系统运行稳定,同化效果良好,计算性能优越,为资料同化系统在并行计算过程中的计算效率提升提供了一种全新的思路。本文设计的优化策略可以保证在不改变同化结果的前提下,使同化性能提高约1倍,大大节约计算时间。此外,这一负载均衡优化方案不仅对于其他使用局地化策略的资料同化方法例如局地集合变换最优插值(Local Ensemble Transform Optimal Interpolation,LETOI)也适用,对大气等其他领域资料同化系统的设计开发也有借鉴意义。
2 算法分析优化及试验方案
2.1 问题介绍与初步优化
首先,对于海洋资料同化而言,陆地网格点由于没有模式数据则无需参与运算。但是,通常的并行计算方案会将所有网格点的计算按CPU核数进行平均划分,即每个核计算网格点的数量相同。这种方案的优势是可以使整个系统简洁、对等和统一,便于程序编写,且更加通用,但会将陆地网格点也纳入其中,影响计算性能。我们以通用地球系统模式(Community Earth System Model,CESM)[16]中约为1°×1°的384×320网格和0.25°×0.25°最优插值海表温度资料(Optimum Interpolation Sea Surface Temperature,OISST)(
图1
图1
模式点与观测点分布图
Fig. 1
Distribution of model grid points and observations
(a) 模式点分布与通常的并行分配方案;(b) OISST观测点分布;经向纬向各只画出1/10的点
(a) Distribution of model grid points and common parallel scheme; (b) Distribution of OISST; Only 1/10 points plotted in horizontal directions
通常的并行方案会采用类似图1a的网格点分配方法,即每个CPU核需要计算的网格点数量相同。考虑到陆地网格点无需计算,被分配到陆地点多的CPU核计算量小,而被分配到陆地点少的CPU核计算量大。由于并行计算中必须等待所有的核计算完成后,才可以进行下一步操作,因此整个程序的效率取决于最慢的核。如果最慢的核耗时过长,无论最快的核计算多快,最终计算性能仍然较差。对于全球海洋资料同化,若依然使用通常的并行方案,则会导致严重的性能问题。
因此,将陆地网格点去除后再将剩余网格点平均分配到每个核上进行计算,是一个较为简单且易实现的方案。我们使用前述的CESM模式及网格,平均选取OISST 5%的观测数据,取局地化半径为3个网格点,使用200核的并行计算环境进行一次LETKF同化,以测试该方案的效果。
图2
2.2 LETKF算法及其时间复杂度分析
LETKF是基于EnKF发展起来的,在EnKF中,若要计算卡尔曼增益矩阵
式中:a,b和o分别表示分析场、背景场、观测,
在实际计算中,为了尽可能减少不必要的计算,通常使用公式变形、存储中间变量等方法,采用如下的计算流程[7]:
(1)对模式变量应用观测算子,计算观测空间扰动
(2)使用集合成员的模式变量计算背景场扰动
(3)对每个网格点,执行局地化操作,选出在网格点局地化半径内的观测数据,同时得到局地空间观测误差协方差
(4)对每个网格点,计算并保存中间变量——局地观测矩阵
(5)对每个网格点,计算集合空间分析场误差协方差的逆
(6)对每个网格点,通过对
(7)计算集合空间卡尔曼增益分析向量
(8)应用卡尔曼增益分析向量,将背景场扰动
我们对该算法的所有步骤进行算法时间复杂度的分析,如表1所示。每一步的算法时间复杂度的具体分析过程,在文中略去。我们仅以第5步为例进行说明。由计算步骤可知,第5步是通过对每个网格点计算
表1 算法时间复杂度分析及说明
Table 1
时间复杂度 | 算法步骤及说明 |
---|---|
(1)应用观测算子,计算观测空间扰动 | |
(2)计算背景场扰动 | |
(3)选择格点内观测数据,得到局地观测误差协方差 | |
(4)计算局地观测中间矩阵 | |
(5)计算集合空间分析场误差协方差的逆 | |
(6)计算集合空间分析场误差协方差与背景场—分析场扰动变换矩阵 | |
(7)计算集合空间卡尔曼增益分析向量 | |
(8)应用卡尔曼增益,得到分析场 |
2.3 LETKF算法瓶颈
通过表1可知,该算法第(5)步与第(6)步的时间复杂度最高,均为4次项,整体算法时间复杂度为
基于以上讨论,我们认为第(5)步是优化的关键,事实上该步骤正是导致计算量分配不均的核心。具体原因在于该步骤执行的是在每个模式网格点的局地化半径内的观测点的相关计算,而模式点或观测点的不均匀分布使得每个模式点局地化半径内的观测点数量出现差异。选择随纬度变化的局地化方案更会加剧这一差异。这一差异在单核计算时不会导致性能问题,但在并行环境下,会导致每个CPU的计算量不相等,使得整体计算性能严重下降。因此,如何尽可能均分计算量是问题的关键。考虑到该步骤中每个网格点的计算量与它局地化半径内的观测点数成正比,则每个网格点的计算量可以用该观测点数来衡量。如果我们能够做到使所有网格的总计算量平均分给每个核进行计算,则可以解决该问题。
将这个问题进行抽象,它等价于将
2.4 使用贪心算法进行优化
贪心算法并不是一种固定的算法,它将一个问题分解成多个连续的小问题,通过解决这一系列小问题来给出大问题的解。其中,在解决每一个小问题时,都做一个对于解决该小问题最佳的决策,期望通过所做的一系列局部最优选择产生出全局的最优解。这种启发式策略常常能给出最优解或非常接近最优的解[25]。
在我们的问题中,大问题为将
根据以上分析,可总结该算法的步骤如下:
(1)计算
(2)将这
(3)对于每个核,做第(4)步操作,得到该核需要计算的网格点集合。
(4)在平衡的二叉搜索树中寻找计算量最接近
由上述步骤可知,该算法是逐核对计算量进行分配,当一个核的分配完成后才进行后面核的分配。在具体的代码实现中,可稍作修改,在每分配一个网格点后,将当前处理的核切换为下一个核,即轮流进行每个核的网格点分配。该修改可使得计算量大的网格先被分配,从而更有可能获得最优解。
需要注意的是,该算法可能会引入额外的开销,需要进行评估。主要有以下2点:
针对第一点,由于该算法是对所有的模式网格点进行一次在二叉搜索树中的查询和修改操作,因此时间复杂度为
因此,应用该算法导致的额外时间消耗几乎可以忽略不计,不会影响到整个同化算法的性能。
2.5 试验方案
基于以上的分析,我们利用贪心算法设计开发了新的基于LETKF的高性能并行资料同化系统,并设计了如下的试验方案,在确保同化的正确性及有效性的前提下,重点关注同化系统的计算性能。
同化系统使用20个集合成员,由模式自由运行1年的结果随机选取,局地化半径在水平方向上为7个网格点,垂直深度为2 000 m,膨胀系数为21%。采用离线方式,独立于模式的并行剖分,在Intel Xeon E5 2.10 GHz处理器集群上使用200核进行同化试验。
同化试验从2004年1月1日开始进行,同化间隔为5天,为期2个月。观测资料采用2004年1月和2月日平均的0.25°×0.25°OISST和中国Argo实时资料中心(www.argo.org.cn)提供的同时期的Argo温度、盐度剖面。所采用的Argo资料的空间分布如图3所示。
图3
图3
2004年1~2月Argo资料空间分布
Fig. 3
Spatial distribution of Argo data in January and February 2004
此外,为了测试贪心算法优化策略的可扩展性,在其他设置不变的情况下,还将使用约为0.1°×0.1°的3 840×3 200高分辨率的网格进行性能测试。
3 试验结果
3.1 结果验证
我们利用同化试验所使用的OISST和Argo资料进行结果验证。
图4
图4
全球海洋SST的RMSE(℃)分布
Fig. 4
Distribution of RMSE(℃)of global ocean SST
(a) 同化前;(b) 同化后
(a) Before data assimilation; (b) After data assimilation
图5
图5
同化前后RMSE垂直分布
Fig. 5
Vertical distribution of RMSE before and after data assimilation
(a) 温度(℃);(b) 盐度;实(虚)线为同化前(后)
(a) Temperature(℃); (b) Salinity; Solid (dashed) line denote RMSE before (after) data assimilation
3.2 性能优化分析
为了分析同化性能优化的效果,我们分别绘制出优化前后每次同化过程中CPU计算的时间分布图。如图6所示,横坐标为CPU核的编号,纵坐标为每个CPU核的计算时间。一次同化所需要的时间由计算时间最长的CPU核决定。
图6
图6
优化前后的CPU计算时间分布
Fig. 6
Distribution of CPU calculation time before and after optimization
(a) 1°×1°分辨率;(b) 0.1°×0.1°分辨率
(a) 1°×1° resolution;(b) 0.1°×0.1° resolution
如图6a,原先每次同化需要约4.2 h,优化之后只需约2.2 h。这只是一次同化所节约的时间,如果进行长时间的同化试验,节约的时间将十分可观。而且如果将来使用更多的观测数据,这一优化策略的效果将更加明显。
另外,该策略在高分辨率(0.1°×0.1°)模式下也具有很好的效果。图6b为0.1°×0.1°的高分辨率网格下CPU计算时间分布图。由于高分辨率模式导致的计算量的增大,优化所节约的时间更加可观,在使用200核的情况下,每次同化可节省约9 d的计算时间。
基于以上结果,该优化策略相比于通常的计算分配策略,可以节约近一半的计算时间,说明该负载均衡策略起到了很好的优化作用。该方法可以使CPU几乎不存在空转的情况,充分利用计算资源,几乎达到了在不修改同化算法情况下的最高计算性能。
4 结论与展望
本文通过对LETKF资料同化方法计算流程的完整分析,研究每一步的算法时间复杂度,从而找到该算法的瓶颈为“计算集合空间分析场误差协方差的逆矩阵”这一步。结合并行计算中CPU计算量分布不均的问题,利用贪心算法提出了一个高效的负载均衡优化方案,设计开发了一个基于LETKF的高性能并行海洋资料同化系统。之后我们将全球海洋2004年1~2月日平均OISST和同时期Argo的温度和盐度剖面作为观测资料进行同化试验,分析了同化的效果和计算的性能。结果显示同化系统工作正常,同化效果良好。与此同时,对比结果表明该系统在保证正确性的前提下,能够提升约1倍的性能,具有在保证二进制一致性即不对结果造成任何影响的情况下,最大化利用计算资源的能力。该优化策略扩展至高分辨率模式后,仍然具有相同的优化能力,证明了系统的可扩展性。由于高分辨率模式计算量的增大,优化所节约的时间变得更加可观。
值得指出的是,该负载均衡方案虽然是基于LETKF,但对其他同样存在计算量分配不均的同化方法(例如LETOI)也可起到很好的优化效果。另外,文中所述的同化系统也可应用于其他海洋(大气)模式,对于大规模的同化试验以及业务化运行的同化系统来说,计算资源消耗和计算成本可以得到显著的降低。该同化系统具有很强的扩展性和移植性,在将来的业务化预报中具有广阔的应用前景。
本文主要以保证二进制一致性为前提,针对LETKF计算量分配不均的问题进行研究与优化。而能否在允许对结果影响较小的前提下针对计算热点本身进行优化还有待进一步研究。
参考文献
The theoretical framework of the ensemble-based data assimilation method and its prospect in oceanic data assimilation
[J].
集合资料同化方法的理论框架及其在海洋资料同化的研究展望
[J]. ,
Current status of global ocean reanalysis datasets
[J].
全球海洋再分析产品的研究现状
[J]. ,
Variational algorithms for analysis and assimilation of meteorological observations: Theoretical aspects
[J]. ,
Development of data assimilation and its application in ocean science
[J].
资料同化技术的发展及其在海洋科学中的应用
[J]. ,
A new approach to linear filtering and prediction problems
[J]. ,
Sequential data assimilation with a nonlinear quasi‐geostrophic model using Monte Carlo methods to forecast error statistics
[J]. ,
Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter
[J]. ,
Applying a four-dimensional local ensemble transform Kalman filter (4D-LETKF) to the JMA nonhydrostatic model (NHM)
[J]. ,
State analysis using the Local Ensemble Transform Kalman Filter (LETKF) and the three-layer circulation structure of the Luzon Strait and the South China Sea
[J]. ,
The Tropical Ocean‐Global Atmosphere observing system: A decade of progress
[J]. ,
Fifteen years of ocean observations with the global Argo array
[J]. ,
Localizing the error covariance by physical distances within a local ensemble transform Kalman filter (LETKF)
[J]. ,
The ensemble Kalman filter in an operational regional NWP system: Preliminary results with real observations
[J]. ,
Ensemble Kalman filter and 4D-Var intercomparison with the Japanese operational global analysis and prediction system
[J]. ,
Introduction to the Community Earth System Model and application to high performance computing
[J].
地球系统模式CESM及其在高性能计算机上的配置应用实例
[J]. ,
Daily high-resolution-blended analyses for sea surface temperature
[J]. ,
The local ensemble transform Kalman filter and the running-in-place algorithm applied to a global ocean general circulation model
[J]. ,
Hindcasts and forecasts of Loop Current and eddies in the Gulf of Mexico using local ensemble transform Kalman filter and optimum-interpolation assimilation schemes
[J]. ,
Data Assimilation of the Global Ocean Using the 4D Local Ensemble Transform Kalman Filter (4D-LETKF) and the Modular Ocean Model (MOM2)
[D].
Solving K-Set Partition Problem Using Genetic Algorithm
[R].
The Differencing Method of Set Partitioning
[M].
Introduction to Algorithms (original book third edition)
[J].
算法导论(原书第3版)
[J]. ,
Diurnal to decadal global forcing for ocean and sea-ice models: The data sets and flux climatologies
[J]. ,
/
〈 | 〉 |