地球科学进展, 2019, 34(5): 531-539 DOI: 10.11867/j.issn.1001-8166.2019.05.0531

研究论文

基于局地集合变换卡尔曼滤波的全球海洋资料同化系统设计及算法加速

范峥,1, 李宏1, 刘向文2, 徐芳华,1

1. 清华大学地球系统科学系,地球系统数值模拟教育部重点实验室,北京 100084

2. 中国气象局国家气候中心,北京,100081

Global Ocean Data Assimilation System Design and Algorithm Acceleration Based on Local Ensemble Transform Kalman Filter

Fan Zheng,1, Li Hong1, Liu Xiangwen2, Xu Fanghua,1

1. Ministry of Education Key Laboratory for Earth System Modeling, Department of Earth System Science, Tsinghua University, Beijing 100084, China

2. National Climate Center, China Meteorological Administration, Beijing 100081, China

通讯作者: 徐芳华(1979-),女,天津人,副教授,主要从事海洋数值模式、数据同化、海洋后报和预报. E-mail:fxu@tsinghua.edu.cn

收稿日期: 2019-01-21   修回日期: 2019-04-02   网络出版日期: 2019-06-11

基金资助: 国家重点研发计划项目“基于高分辨率气候系统模式的无缝隙气候预测系统研制与评估” .  编号:2016YFA0602100

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

FanZheng(1994-),male,BengbuCity,AnhuiProvince,Masterstudent.Researchareasincludehighperformancecomputinginearthscience.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倍,说明新设计的负载均衡方案稳定可靠。该方案具有很强的可扩展性和移植性,在业务预报中有广泛的应用前景。

关键词: 海洋资料同化 ; 局地集合变换卡尔曼滤波 ; 负载均衡策略 ; 高性能

Abstract

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: Ocean data assimilation ; Local Ensemble Transform Kalman filter ; Load balancing strategy ; High performance

PDF (2535KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

范峥, 李宏, 刘向文, 徐芳华. 基于局地集合变换卡尔曼滤波的全球海洋资料同化系统设计及算法加速. 地球科学进展[J], 2019, 34(5): 531-539 DOI:10.11867/j.issn.1001-8166.2019.05.0531

Fan Zheng, Li Hong, Liu Xiangwen, Xu Fanghua. Global Ocean Data Assimilation System Design and Algorithm Acceleration Based on Local Ensemble Transform Kalman Filter. Advances in Earth Science[J], 2019, 34(5): 531-539 DOI:10.11867/j.issn.1001-8166.2019.05.0531

1 引 言

资料同化是数值天气预报的关键技术之一。它不仅可以有效地将时空分布不均匀、不同来源及不同观测误差的常规和非常规观测大气—海洋资料有效地融入数值模式[1],获得更优的状态变量,从而能够为数值预报模式提供更准确更丰富的初值,同时也可以制作高时空分辨率的再分析资料集,还能为数值模式相关物理量和参数的设置提供依据[2]。因此,资料同化在大气—海洋科学的研究中占有非常重要的地位。

然而,高昂的计算代价始终是制约资料同化技术发展的关键因素,资料同化方法始终在平衡同化效果和计算量之间发展。这主要体现在同化方法和观测资料两方面。在同化方法方面,例如连续资料同化方法中的三维变分(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)(https://www.ncdc.noaa.gov/oisst)[17]为例,画出模式点与观测点分布图(图1)。

图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同化,以测试该方案的效果。

去掉陆地点前后CPU计算时间分布如图2所示。从图2圆形点可以看出,正是由于陆地点的存在导致CPU计算量分配不均,使得不同核之间的计算时间差距过大。如图2所示,在去掉陆地点前,最慢的核计算耗时约150 s,而在去掉陆地后,最慢的核计算耗时约96 s,节省了约36%的计算时间。

图2

图2   CPU计算时间分布

Fig. 2   Distribution of CPU calculation time


虽然这种方案已经使得性能得到了提升,但是观察图2的三角形点不难发现,最快与最慢的核之间仍然相差约45 s,这说明不同核之间的计算量仍然存在差距。如果我们能够消除这一差距,使图2中CPU计算时间的分布更接近于一条直线,则会进一步提升性能。在大规模观测数据参与同化计算的情况下,这一优化至关重要,每次同化都可节约数小时甚至数天的时间。欲实现这一目标,需要对LETKF算法进行详尽的分析,分析该算法每个步骤的时间复杂度,找到成为瓶颈的关键步骤,分析其中CPU计算量不均等的原因并加以解决。

2.2 LETKF算法及其时间复杂度分析

LETKF是基于EnKF发展起来的,在EnKF中,若要计算卡尔曼增益矩阵K,需要在观测空间上计算逆矩阵,这个计算量通常是不可接受的。相比于EnKF,LETKF采用局地化策略即选取一个局地化半径,只针对该半径之内的观测数据进行同化,同时运用集合变换,从而在每个点的集合空间上计算其逆矩阵。LETKF通过以下方程的求解得到分析场[7]

Yb=i=1kH(xb(i))-H(xb)¯P˜a=(Yb)TR-1Yb+(k-1)I/ρ-1K=XbP˜a(Yb)TR-1x¯a=x¯b+Kyo-H(xb)¯Xa=Xb(k-1)P˜a1/2xa=x¯a+Xa

式中:a,b和o分别表示分析场、背景场、观测,Yb是观测空间扰动,Xb是背景场扰动,R是局地空间观测误差协方差矩阵,P˜a是集合空间分析场误差协方差,K是卡尔曼增益矩阵,H是观测算子,k是集合数。

在实际计算中,为了尽可能减少不必要的计算,通常使用公式变形、存储中间变量等方法,采用如下的计算流程[7]

(1)对模式变量应用观测算子,计算观测空间扰动YbYbk个列向量yb(i)-y¯b组成,yx执行观测算子得到。

(2)使用集合成员的模式变量计算背景场扰动XbXbk个列向量xb(i)-x¯b组成。

(3)对每个网格点,执行局地化操作,选出在网格点局地化半径内的观测数据,同时得到局地空间观测误差协方差R

(4)对每个网格点,计算并保存中间变量——局地观测矩阵C=(Yb)TR-1

(5)对每个网格点,计算集合空间分析场误差协方差的逆(P˜a)-1=CYb+(k-1)I/ρ。该步骤并非求逆,而是将集合空间分析场误差协方差的逆作为一个整体。

(6)对每个网格点,通过对(P˜a)-1求逆矩阵得到集合空间分析场误差协方差P˜a,并计算背景场—分析场扰动变换矩阵Wa=[(k-1)P˜a]1/2,以变换至集合空间。

(7)计算集合空间卡尔曼增益分析向量{wa(i)},通过均值w¯a=P˜aC(yo-y¯b)与扰动变换矩阵Wa求和得到。

(8)应用卡尔曼增益分析向量,将背景场扰动Xbwa(i)相乘后与x¯b求和,得到分析场{xa(i)}

我们对该算法的所有步骤进行算法时间复杂度的分析,如表1所示。每一步的算法时间复杂度的具体分析过程,在文中略去。我们仅以第5步为例进行说明。由计算步骤可知,第5步是通过对每个网格点计算CYb+(k-1)I/ρ得到集合空间分析场误差协方差的逆(P˜a)-1。其中,C是一个k×li的矩阵,Yb是一个li×k的矩阵,CYb矩阵相乘将得到一个k×k的矩阵。该矩阵的每个元素的计算需要li次相乘求和操作。因此,计算CYb需要k2li次操作,时间复杂度为Ok2li。计算(k-1)I/ρ并与CYb求和需要k2次操作,时间复杂度为Ok2。总体时间复杂度为Ok2li+k2=Ok2li。由于所有m[g]个网格点均需该操作,可得该步骤的时间复杂度则为Ok2i=1m[g]li

表1   算法时间复杂度分析及说明

Table 1  Algorithm time complexity analysis and

时间复杂度算法步骤及说明
Okl[g](1)应用观测算子,计算观测空间扰动
O(km[g])(2)计算背景场扰动
Oki=1m[g]li(3)选择格点内观测数据,得到局地观测误差协方差
Oki=1m[g]li(4)计算局地观测中间矩阵
Ok2i=1m[g]li(5)计算集合空间分析场误差协方差的逆
Ok3m[g](6)计算集合空间分析场误差协方差与背景场—分析场扰动变换矩阵
Oki=1m[g]li+k(7)计算集合空间卡尔曼增益分析向量
Ok2m[g](8)应用卡尔曼增益,得到分析场

注:k为集合数,m[g]为模式网格点数,l[g]为观测点数,li为每个模式网格点的局地空间观测点数

新窗口打开| 下载CSV


2.3 LETKF算法瓶颈

通过表1可知,该算法第(5)步与第(6)步的时间复杂度最高,均为4次项,整体算法时间复杂度为Ok2i=1m[g]li+k3m[g],瓶颈在哪一步取决于lik的关系。对海洋资料同化而言,通常使用20~40个集合成员即可取得较好的效果[9,18,19]li的大小取决于局地化半径的选取,一般选择固定的5~7个网格点的距离[10,19]或根据纬度在2°~10°线性变化[20]。以选定20个集合成员、局地化半径为5个网格点、使用0.25°×0.25°的OISST观测资料为例,此时li=(50.25×2)2=2 500,远大于k=20,整体算法时间复杂度为Ok2i=1m[g]li+k3m[g]=Ok2i=1m[g]li,整个算法的瓶颈在于第(5)步。在实际情况中也可能存在只有部分网格点lik成立的情况,但只要在并行计算时存在某个核的l¯k,第(5)步即成为瓶颈。若在某些应用中既无如此大量的观测数据,又不存在一个观测数据的热点区域,则此时整个算法本身的计算开销已经很低,提升性能的必要性大大降低。

基于以上讨论,我们认为第(5)步是优化的关键,事实上该步骤正是导致计算量分配不均的核心。具体原因在于该步骤执行的是在每个模式网格点的局地化半径内的观测点的相关计算,而模式点或观测点的不均匀分布使得每个模式点局地化半径内的观测点数量出现差异。选择随纬度变化的局地化方案更会加剧这一差异。这一差异在单核计算时不会导致性能问题,但在并行环境下,会导致每个CPU的计算量不相等,使得整体计算性能严重下降。因此,如何尽可能均分计算量是问题的关键。考虑到该步骤中每个网格点的计算量与它局地化半径内的观测点数成正比,则每个网格点的计算量可以用该观测点数来衡量。如果我们能够做到使所有网格的总计算量平均分给每个核进行计算,则可以解决该问题。

将这个问题进行抽象,它等价于将n个整数{a1,,an}分成k组,要求每组数之和相等。在我们的问题中,n就是模式网格点数,ai为每个网格点的观测空间大小,k为使用的CPU核数。该问题是计算机领域经典的K划分问题(K-Partition Problem),它属于计算复杂性理论中的非确定性多项式问题(Non-deterministic Polynomial-time Hard,NP-Hard),在理论上已经证明不存在时间复杂度为多项式的解法,需要至少指数级别的时间[21],无法在可以接受的时间内得到精确解。然而在我们的应用中,我们并不需要得到该问题的精确解,只需求得一个相对较好的近似解。

2.4 使用贪心算法进行优化

目前,该问题有很多近似解法,例如贪心算法[22]、遗传算法[23]和最大差分法[24]等。相比于遗传算法和最大差分法,贪心算法最为简单、易于实现,且效果较好,因此采用该算法来求解。

贪心算法并不是一种固定的算法,它将一个问题分解成多个连续的小问题,通过解决这一系列小问题来给出大问题的解。其中,在解决每一个小问题时,都做一个对于解决该小问题最佳的决策,期望通过所做的一系列局部最优选择产生出全局的最优解。这种启发式策略常常能给出最优解或非常接近最优的解[25]

在我们的问题中,大问题为将n个整数{a1,,an}分成k组,要求每组数之和近似相等,其值约为S¯=1ki=1nai。我们首先可将该问题拆分成k个子问题,即对于其中某一组,如何从当前剩下的数中取出一些数,使得该组数之和尽可能接近S¯。接着,对该子问题进一步拆分,即若此时该组数之和为S0,如何从当前剩下的数中取出一个数使得该数尽可能接近S¯-S0

我们分别将这2个子问题命名为问题1和问题2。通过多次解决问题2即可解决问题1;通过k次解决问题1即可解决原问题。因此,问题的关键在于如何解决问题2。将问题2换一种表述,它等价于如何从m个整数{a1,,am}中找出一个最接近于S的数并删除。对于这一问题,解决方法可以非常简单,只需进行一次遍历,其时间复杂度为Om。但是,为了不造成额外的开销,必须尽可能使用最快的算法。为了解决这一关键问题,需要借助一些特殊的数据结构,通常使用平衡的二叉搜索树[25]。该数据结构的特点是可以在O(logm)的时间内进行查询和修改操作,使得时间复杂度大大降低。平衡的二叉搜索树在具体实现上主要有2种常见方式:AVL-Tree和RB-Tree[25]

根据以上分析,可总结该算法的步骤如下:

(1)计算S¯=1ki=1nai,即期望的在最好情况下每个核的计算量。

(2)将这n个网格点的计算量ai组织成一个平衡的二叉搜索树,相当于n次插入操作。

(3)对于每个核,做第(4)步操作,得到该核需要计算的网格点集合。

(4)在平衡的二叉搜索树中寻找计算量最接近S¯-si的网格,将其加入该核的网格集合并在二叉搜索树中删除,直到siS¯即该核心已达到期望被分配的计算量。

由上述步骤可知,该算法是逐核对计算量进行分配,当一个核的分配完成后才进行后面核的分配。在具体的代码实现中,可稍作修改,在每分配一个网格点后,将当前处理的核切换为下一个核,即轮流进行每个核的网格点分配。该修改可使得计算量大的网格先被分配,从而更有可能获得最优解。

需要注意的是,该算法可能会引入额外的开销,需要进行评估。主要有以下2点:贪心算法本身的计算开销;修改并行剖分导致的数据传输开销。

针对第一点,由于该算法是对所有的模式网格点进行一次在二叉搜索树中的查询和修改操作,因此时间复杂度为O(nlogn)n为模式网格点的个数。由于目前的CPU的单核性能可以达到每秒数十亿次计算,因此贪心算法几乎没有执行开销;针对第二点,需要对传输的数据量和带宽进行分析。该算法要求重新传输所有集合的所有网格的数据,数据量约为集合数×模式网格点数×垂向层数×变量所占字节数。以温度为例,数据量约为20×320×384×60×4590 MB。计算集群使用四路InfiniBand FDR,其带宽为56 Gb/s,即每秒传输7 168 MB的数据,数据量仅占很小的带宽。

因此,应用该算法导致的额外时间消耗几乎可以忽略不计,不会影响到整个同化算法的性能。

2.5 试验方案

基于以上的分析,我们利用贪心算法设计开发了新的基于LETKF的高性能并行资料同化系统,并设计了如下的试验方案,在确保同化的正确性及有效性的前提下,重点关注同化系统的计算性能。

同化试验所采用的模式是CESM中的并行海洋模块2(The Parallel Ocean Program 2,POP2)全球海洋模式。水平网格约为1°×1°,在赤道附近加密至约0.25°,至两极逐渐减小至1.875°(图1a)。水平格点数为384×320,垂向分为60层。在海气界面间,海洋模式受到COREv2提供的大气气候态风场、海面压强、气温、蒸发和降水等驱动[26]

同化系统使用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   20041~2Argo资料空间分布

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资料进行结果验证。

4a和4b分别为同化前后全球海洋海表温度(Sea Surface Temperature,SST)均方根误差(Root Mean Square Error, RMSE)的空间分布图,图5则给出同化前后温度和盐度的RMSE垂直分布。

图4

图4   全球海洋SSTRMSE)分布

Fig. 4   Distribution of RMSEof 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


结果表明,同化前,全球海洋SST的RMSE普遍大于1 ℃(图4a),同化后,SST的RMSE普遍小于0.7 ℃(图4b)。同化前后的温度和盐度RMSE廓线表明,同化前,温度(盐度)廓线的RMSE最大可达1.9 ℃(0.75),同化后,对应最大RMSE则降为0.8 ℃(0.4)以下。同化有效修正了模式的温度和盐度结果,同化效果明显。

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计算量分配不均的问题进行研究与优化。而能否在允许对结果影响较小的前提下针对计算热点本身进行优化还有待进一步研究。

参考文献

Shen Zheqi , Tang Youmin , Gao Yanqiu .

The theoretical framework of the ensemble-based data assimilation method and its prospect in oceanic data assimilation

[J]. Acta Oceanologica Sinica, 2016, 38(3):1-14.

[本文引用: 1]

沈浙奇, 唐佑民, 高艳秋 .

集合资料同化方法的理论框架及其在海洋资料同化的研究展望

[J]. 海洋学报, 2016, 38(3):1-14.

[本文引用: 1]

Wang Shihong , Zhao Yiding , Yin Xunqiang , et al .

Current status of global ocean reanalysis datasets

[J]. Advances in Earth Science, 2018, 33(8):794-807.

[本文引用: 1]

王世红,赵一丁,尹训强, .

全球海洋再分析产品的研究现状

[J]. 地球科学进展, 2018, 33(8):794-807.

[本文引用: 1]

Le Dimet F X , Talagrand O .

Variational algorithms for analysis and assimilation of meteorological observations: Theoretical aspects

[J]. Tellus A: Dynamic Meteorology and Oceanography, 1986, 38(2):97-110.

[本文引用: 1]

Li Hong , Xu Jianping .

Development of data assimilation and its application in ocean science

[J]. Marine Science Bulletin, 2011, 30(4):463-472.

[本文引用: 1]

李宏, 许建平 .

资料同化技术的发展及其在海洋科学中的应用

[J]. 海洋通报, 2011, 30(4):463-472.

[本文引用: 1]

Kalman R E .

A new approach to linear filtering and prediction problems

[J]. Journal of Fluids Engineering, 1960811):35-45.

[本文引用: 1]

Evensen G .

Sequential data assimilation with a nonlinear quasi‐geostrophic model using Monte Carlo methods to forecast error statistics

[J]. Journal of Geophysical Research: Oceans, 1994, 99(C5):10 143-10 162.

[本文引用: 1]

Hunt B R , Kostelich E J , Szunyogh I .

Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter

[J]. Physica D, 2007, 230(1/2):112-126.

[本文引用: 3]

Miyoshi T , Aranami K .

Applying a four-dimensional local ensemble transform Kalman filter (4D-LETKF) to the JMA nonhydrostatic model (NHM)

[J]. Sola, 2006, 2: 128-131.

[本文引用: 1]

Miyoshi T .Ensemble data assimilation for idealized California current system with ROMS-LETKF[C]//14th Symposium on Integrated Observing and Assimilation Systems for the Atmosphere, Oceans, and Surface Land . 2010.

[本文引用: 2]

Xu F H , Oey L Y .

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]. Ocean Dynamics, 2014, 64(6):905-923.

[本文引用: 2]

Mcphaden M J , Busalacchi A J , Cheney R , et al .

The Tropical Ocean‐Global Atmosphere observing system: A decade of progress

[J]. Journal of Geophysical Research: Oceans, 1998,103(C7). DOI:10.1029/97JC0926 .

[本文引用: 1]

Riser S C , Freeland H J , Roemmich D , et al .

Fifteen years of ocean observations with the global Argo array

[J]. Nature Climate Change, 2016, 6(2):145-153.

[本文引用: 1]

Miyoshi T , Yamane S , Enomoto T .

Localizing the error covariance by physical distances within a local ensemble transform Kalman filter (LETKF)

[J]. Sola, 2007, 3: 89-92.

[本文引用: 1]

Bonavita M , Torrisi L , Marcucci F .

The ensemble Kalman filter in an operational regional NWP system: Preliminary results with real observations

[J]. Quarterly Journal of the Royal Meteorological Society, 2008, 134(636): 1 733-1 744.

[本文引用: 1]

Miyoshi T , Sato Y , Kadowaki T .

Ensemble Kalman filter and 4D-Var intercomparison with the Japanese operational global analysis and prediction system

[J]. Monthly Weather Review, 2010, 138(7): 2 846-2 866.

[本文引用: 1]

Wan Xiuquan , Liu Zedong , Shen Biao , et al .

Introduction to the Community Earth System Model and application to high performance computing

[J]. Advances in Earth Science, 2014, 29(4):482-491.

[本文引用: 1]

万修全,刘泽栋,沈飙, .

地球系统模式CESM及其在高性能计算机上的配置应用实例

[J]. 地球科学进展, 2014, 29(4):482-491.

[本文引用: 1]

Reynolds R W , Smith T M , Liu C , et al .

Daily high-resolution-blended analyses for sea surface temperature

[J]. Journal of Climate, 2007, 20(22):5 473-5 496.

[本文引用: 1]

Penny S G , Kalnay E , Carton J A , et al .

The local ensemble transform Kalman filter and the running-in-place algorithm applied to a global ocean general circulation model

[J]. Nonlinear Processes in Geophysics, 2013, 20(6):1 031-1 046.

[本文引用: 1]

Xu F H , Oey L Y , Miyazawa Y , et al .

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]. Ocean Modelling, 2013, 69: 22-38.

[本文引用: 2]

Penny S G .

Data Assimilation of the Global Ocean Using the 4D Local Ensemble Transform Kalman Filter (4D-LETKF) and the Modular Ocean Model (MOM2)

[D]. Maryland:University of Maryland, 2011.

[本文引用: 1]

Hayes B .

Computing science: The easiest hard problem

[J]. American Scientist, 2002, 90(2):113-117.

[本文引用: 1]

Korf R E .

Multi-way number partitioning

[C]//International Joint Conference on Ijcai. DBLP, 2009.

[本文引用: 1]

Ahmed M S , Nisha T A .

Solving K-Set Partition Problem Using Genetic Algorithm

[R]. Chicago:East West University, 2014.

[本文引用: 1]

Karmarker N , Karp R M .

The Differencing Method of Set Partitioning

[M]. Campbell Hall:Berkeley California University of California at Berkeley, 1983.

[本文引用: 1]

Cormen T H , Leiserson C E , Rivest R L , et al .

Introduction to Algorithms (original book third edition)

[J]. Computer Education, 2013(10):237, 161-192.

[本文引用: 3]

Cormen T H , Leiserson C E , Rivest R L , .

算法导论(原书第3版)

[J]. 计算机教育, 2013(10):237, 161-192.

[本文引用: 3]

Large W G , Yeager S G .

Diurnal to decadal global forcing for ocean and sea-ice models: The data sets and flux climatologies

[J]. National Center for Atmospheric Research, 2004. DOI:10.5065/D6KK98Q6 .

[本文引用: 1]

/