GRAPES全球模式静力平衡奇异向量改进及应用试验*

霍振华 李晓莉 陈 静 刘永柱

1.中国气象局数值预报中心，北京，100081

2.国家气象中心，北京，100081

1 引 言

由于初始误差、模式误差的存在，以及大气运动本身的混沌特性，单一确定性预报存在不确定性，集合预报正是解决这一问题的有效手段（Leith，1974）。作为制作集合预报的关键，初始扰动的生成极为重要，其始于20世纪70年代。几十年来，国际上已发展了多种集合预报初值扰动方法。其中，奇异向量方法（Singular Vector，SV；Lorenz，1965；Mureau， et al， 1993； Molteni， et al， 1996； Buizza，1997；Leutbecher，et al，2008）已成功应用于欧洲中期 天 气 预 报 中 心 （ European Centre for Medium-Range Weather Forecasts，ECMWF）、日本气象厅（Japan Meteorological Agency，JMA）、加拿大气象局（Canadian Meteological Center，CMC）和中国气象局（China Meteological Center，CMA）。

奇异向量（SV）基于非线性模式的切线性和伴随模式计算得到，是线性模式中增长最快的初始扰动，可以描述大气在中高纬度的斜压不稳定（Buizza，et al，1993；Hoskins，et al，2000）。此外，奇异向量方法在热带气旋（Tropical cyclone，TC）集合预报方面具有较好的表现。热带地区的非绝热过程很重要且具有很强的非线性，对于热带地区，奇异向量对切线性模式中非绝热物理过程很敏感，如果计算整个热带地区的奇异向量（简称热带SV），有时会生成在非线性模式中不增长的虚假结构（Barkmeijer，et al，2001）。因此，对于热带地区，当存在热带气旋时，ECMWF全球集合预报系统计算针对热带气旋目标区的湿（切线性和伴随模式中使用非绝热物理过程）奇异向量（即湿TCSV），以表征热带气旋内部及附近区域的初始不确定性（Barkmeijer，et al，2001；Puri，et al，2001），从而对热带气旋预报比较有效。Zadra等（2004）指出，奇异向量对线性物理过程非常敏感，使用不同的线性物理过程组合可以得到具有不同结构和发展特征的奇异向量。Kim等（2009）指出，湿线性物理过程和湿总能量范数可以增大奇异向量的增长率，使其具有更小尺度的水平结构和更加集中在边界层的扰动分布，使用湿线性物理过程和湿总能量范数得到的奇异向量其总能量大值集中在热带气旋周围。研究不同物理过程对奇异向量空间结构及其增长特征的影响，不仅可以为热带气旋集合预报初值扰动的生成提供科学依据，还可以帮助分析线性物理过程的合理性。

全球区域一体化同化预报系统（Global/Regional Assimilation and Prediction System， GRAPES）是在中华人民共和国科学技术部和中国气象局支持下中国自主研发的新一代数值预报系统（陈德辉等，2006；薛纪善等，2008；Zhang，et al，2008；龚建东等，2016；苏勇等，2018；刘艳等，2019）。基于GRAPES 全球预报系统（GRAPES Global Forecast System，GRAPES-GFS）的切线性和伴随模式，中国气象局数值预报中心开展了基于总能量模的全球奇异向量求解技术及基于奇异向量产生集合预报初值扰动的研发（刘永柱等，2011，2013；李晓莉等，2019a，2019b；霍振华等，2020），并在 2018年 12月业务运行的GRAPES全球集合预报系统得到应用。如前所述，切线性和伴随模式中湿线性化物理过程方案的使用对TCSV的结构有重要影响，目前在GRAPES_GEPS中，与中高纬度目标区奇异向量计算设置相同，TCSV计算也采用干线性物理过程的设置（线性边界层（Planetary Boundary Layer，PBL）方案和线性次网格尺度地形拖曳（Mountain Blocking，MB）方案），这会造成TCSV的结构在某种程度上不能较好地体现热带气旋初值的不确定性。随着GRAPES-GFS切线性和伴随模式的完善，包括了更多的线性化物理过程（刘永柱等，2019），尤其是湿线性物理过程如大尺度凝结和积云对流，这使得基于GRAPES-GFS开展不同线性物理过程对奇异向量，尤其是TCSV的影响研究成为可能。然而，在研究这些线性化物理过程对TCSV的影响时，发现存在部分奇异向量非线性积分崩溃的问题。这些非线性积分崩溃的奇异向量，一方面说明其结构可能存在扰动变量之间不协调之处，另一方面阻碍分析线性物理过程对奇异向量的影响，从而影响到奇异向量初值扰动方法的优化改进。因此，有必要对原有的奇异向量求解方法进行改进，以解决奇异向量非线性积分崩溃问题，进而改进GRAPES全球集合预报效果。

GRAPES全球数值预报模式（薛纪善等，2008）是一个非静力、半隐式半拉格朗日的均匀经纬度格点模式，采用了高度地形追随坐标。为了提高垂直运动方程的计算精度，模式引入含复杂地形的参考大气，参考大气温度为常数（T0），参考Exner气压变量（Πr）和参考位温（θr）满足静力平衡关系。GRAPESGFS 在近些年得到了快速的发展（Huang Y Y，et al，2013；Huang B，et al，2014），可以较好地预报副热带高压（Chen，et al，2008）。

在GRAPES全球数值预报模式的基础上，中国气象局数值预报中心发展了预报时效为15 d的GRAPES全球集合预报系统（GRAPES-GEPS），集合成员数为 31，水平分辨率为 0.5°×0.5°，垂直层数为60。GRAPES-GEPS基于奇异向量方法产生初值扰动（刘永柱等，2011，2013；李晓莉等，2019a，2019b；霍振华等，2020），使用随机物理倾向项扰动方 法 （Stochastically perturbed parameterization tendencies，SPPT；袁月等，2016；李晓莉等，2019a）和动能后向散射补偿扰动方法（Stochastic Kinetic Energy Backscatter，SKEB；彭飞等，2019）对模式物理过程倾向进行随机扰动。作为GRAPES-GEPS的核心，奇异向量求解采用分目标区计算，分别以北半球副热带地区、南半球副热带地区为目标区域计算奇异向量，进而对其进行线性组合以产生初值扰动，并将初值扰动以正、负扰动对的形式叠加在分析场上生成扰动初值场。当存在热带气旋时，以热带气旋为中心10°×10°的范围为目标区域求解TCSV，将TCSV和副热带SV及演化SV线性组合产生初值扰动。

针对TCSV在非线性积分时出现计算崩溃问题，为研究线性物理过程对奇异向量的影响创造条件，有必要改进奇异向量求解方法，进而改善GRAPES全球集合预报效果。基于已经建立的GRAPES-GEPS，改进奇异向量求解方法，发展了静力平衡奇异向量改进方法，并进行集合预报试验，对比分析奇异向量求解方法改进前、后TCSV的结构特征和初值扰动特征，评估静力平衡奇异向量改进方法对北半球和南半球等压面要素集合预报、热带气旋路径集合预报和中国地区24 h累计降水概率预报技巧的影响。

2 试验方案与资料

2.1 试验方案

2.1.1 原有的奇异向量求解方法

GRAPES全球模式奇异向量的计算需要使用相应的切线性和伴随模式。在计算奇异向量时，忽略初始时刻的垂直运动，GRAPES切线性模式的大气基本状态向量为（u，v，θ′，Π′），其扰动向量为GRAPES模式预报变量的扰动向量，分别由纬向风扰动向量（u′）、经向风扰动向量（v′）、扰动位温的扰动向量 (θ′)′和扰动Exner气压变量的扰动向量(Π′)′组成，即扰动向量x=(u′,v′,(θ′)′,(Π′)′)。有关奇异向量求解的细节参见刘永柱等（2011）。下面简要介绍一下奇异向量的计算方案。

对于给定的权重模，奇异向量是在最优化时间间隔内增长最快的一组正交扰动。对于一个初始扰动向量x0，经过GRAPES全球切线性模式（记为L）一定时间的向前积分可以获得线性演化的扰动向量xt=Lx0。在本文中，取最优化时间间隔为2 d，求解奇异向量和演化SV。GRAPES全球奇异向量的求解可归结为演化扰动向量(xt)模与初始扰动向量(x0)模比值的最大化问题。

式中，λ为奇异向量值，E为度量扰动大小的权重模，P是投影算子，上标T是转置符号。为获得某个特定目标区域内增长最快的奇异向量扰动，在此引入了投影算子P，通过该投影算子将实际计算过程中目标区之外的扰动设置为0。式（1）可转换为如下的奇异值分解问题。

式中，xi表示求解出的第i个奇异向量， λi为对应的第i个奇异值。式中，LT为GRAPES全球切线性模式的转置，也就是伴随模式。对于式（2）定义的GRAPES奇异向量计算问题可以采用洛伦兹迭代算法求解，在迭代过程中需多次积分切线性模式L及伴随模式LT。

在计算GRAPES奇异向量时，采用的权重模为总能量模（刘永柱等，2013），扰动量的总能量模（E）的计算如下

式中，前面两项之和为扰动动能（kinetic energy，KE）模，后两项之和表示扰动位能（potential energy，PE）模，其中第3项和第4项分别表示位能中扰动位 温 (θ′)′和 扰 动 无 量 纲 气 压 (Π′)′分 量 的 贡 献 。 式（3）中，为地形追随坐标，λ和 φ 分别代表模式球面坐标系中的经度和纬度，cp为干空气的定压比热。Tr、θr、Πr和 ρr分别表示参考温度、参考位温、参考无量纲气压和参考密度，其计算公式分别为

式中，Rd为干空气气体常数，pr为标准大气，Tr为参考温度，g为重力加速度。

以2019年8月8日12时（世界时，下同）为初始时刻，此时存在1909号台风“利奇马”和1910号台风“罗莎”。表1给出了不同试验方案的配置，其中PBL为线性边界层方案，MB为线性次网格尺度地形拖曳方案，LC为线性大尺度凝结方案，CU为线性积云对流方案。针对不同的试验方案，以2019年8月8日12时为初始时刻，采用原有的奇异向量计算方法，求解了两个台风对应的前5个TCSV。针对表1中的5组试验方案，分别进行了计算出的奇异向量非线性积分试验（积分48 h），以分析奇异向量扰动，特别是湿线性化物理过程奇异向量在非线性模式积分中的稳定性。结果表明，试验方案test_pbl、test_lc、test_cu、test_lc_cu两个台风都存在对应的TCSV非线性积分崩溃的情况。test_mb试验第21个北半球奇异向量同样出现了非线性积分崩溃的现象。这些非线性积分崩溃的奇异向量，一方面说明其结构可能存在扰动变量之间不协调之处，另一方面阻碍分析线性物理过程对奇异向量的影响，从而影响到奇异向量初值扰动方法的优化改进。因此，需要改进奇异向量求解方法，解决其非线性积分崩溃问题。

表1 不同试验方案所使用的线性物理过程Table 1 Linearized physical processes used for different experiment schemes

2.1.2 静力平衡奇异向量改进方法

原奇异向量求解方法初始时刻，扰动位温(θ′)′和扰动无量纲气压 (Π′)′分量是式（2）中需求解的奇异向量的两个独立扰动变量分量。在洛伦兹迭代算法流程中，独立进行切线性模式和伴随模式循环积分迭代，最终获得独立的奇异向量扰动位温和扰动无量纲气压分量。上述处理可能导致扰动初始场的静力平衡关系（式（5））遭到较为严重破坏，从而使奇异向量在非线性模式积分时崩溃。奇异向量非线性积分崩溃，一方面说明其结构可能存在扰动变量之间不协调之处，另一方面阻碍分析线性物理过程对奇异向量的影响，从而影响到奇异向量初值扰动方法的优化改进。为此，本研究对奇异向量计算流程（刘永柱等，2013）进行改进，使得扰动位温和扰动无量纲气压分量之间符合静力平衡关系。在迭代求解的每一步，对初始扰动向量（进行切线性模式积分的扰动向量）中的气压扰动分量根据静力平衡关系，由扰动向量中的位温扰动和比湿扰动来获得（如式（6）所示），其他计算流程和原奇异向量求解方法保持一致，进而循环迭代获得扰动变量具有静力平衡关系特征的GRAPES SV。

式中，Π为无量纲气压，q为比湿，θ为位温。

式中，q′为比湿扰动。目前，在GRAPES全球集合预报业务系统中，初值扰动尚未包括比湿扰动，奇异向量不包括比湿扰动。因此，在本研究中，q′=0。

基于静力平衡奇异向量改进方法，分别重新计算两个台风对应于表1中不同试验方案的前5个TCSV，并进行了对应的非线性模式积分。结果表明，使用静力平衡奇异向量改进方法得到的TCSV，均可以正常非线性积分。此外，在新方案中采用test_mb方案时，北半球和南半球奇异向量都可以正常非线性积分，在原奇异向量计算方案中北半球第21个奇异向量非线性积分崩溃的问题也得到了解决。

静力平衡奇异向量改进方法可以解决奇异向量的非线性积分崩溃问题，然而奇异向量求解方法的最终目的是产生集合预报所需要的初值扰动以期得到具有较高技巧的集合预报结果，因此有必要分析该方法对奇异向量结构、初值扰动及集合预报技巧的影响。目前GRAPES全球集合预报业务系统中使用的线性物理过程为PBL和MB，因此采用线性物理过程PBL和MB对比分析奇异向量求解方法改进前、后TCSV和初值扰动的结构特征及其集合预报技巧。原奇异向量求解方法简称为原方法，静力平衡奇异向量改进方法简称为新方法。

2.1.3 试验设计

为分析新方法对GRAPES-GEPS奇异向量及初值扰动的影响，选取2019年8月8日12时为初始时刻，采用原方法和新方法，分别求解台风利奇马和罗莎对应的TCSV，分析了两种方法对应的TCSV和初值扰动的结构特征。进而开展了基于两种方法初值扰动的集合预报试验（2019年8月4、6、8、11、13和 15日每天进行，12时起报），分析了两组集合预报试验对北半球和南半球等压面要素集合预报、热带气旋路径集合预报和中国地区24 h累计降水概率预报技巧的差异。在进行集合预报试验时，采用SPPT和SKEB模式扰动方案，集合成员数为31个。

2.2 资 料

基于GRAPES-GEPS求解奇异向量并进行集合预报试验时，初值根据GRAPES-GFS业务同化分析场插值得到，初始时刻热带气旋定位信息来自国家气象信息中心GTS通讯网络数据，台风路径和强度集合预报检验评估采用中央气象台提供的热带气旋最佳路径数据，降水预报检验采用国家气象信息中心2400多个国家级自动气象观测站的降水观测资料。

3 静力平衡奇异向量改进方法对奇异向量和初值扰动结构的影响

选取2019年8月8日12时为初始时刻，采用原方法和新方法，分别求解台风利奇马和罗莎对应的TCSV，分析了两种方法对应的TCSV和初值扰动的结构特征。由于台风罗莎和利奇马对应的TCSV结构类似，以下以台风利奇马为例，分析两种方法对应的TCSV结构特征。

3.1 热带气旋奇异向量结构分析

图1给出了原方法和新方法台风利奇马对应的前5个TCSV（放大1000倍）对应的能量分布特征，其中，图1a1—a3和图1b1—b3分别对应的是原方法和新方法对应的TCSV。图1a1、b1对应的是不同TCSV对应的动能和位能分布，其中KE为动能，PE为位能；图1a2、b2对应的是5个TCSV动能和位能及总能量均值的垂直分布，其中TE为总能量；图1a3、b3对应的是5个TCSV动能和位能及总能量的能量谱均值分布。结果表明，无论采用原方法还是新方法求解奇异向量，TCSV具有类似的能量垂直和空间分布特征。

图2给出了原方法和新方法台风利奇马对应的第3个TCSV（放大1000倍）的位温扰动分量在模式第26层（大约500 hPa）的水平结构和在25°N的垂直结构。可以看出，两种奇异向量计算方法获得的TCSV中位温扰动无论是初始时刻和线性演化位温扰动的水平分布和垂直结构是非常相似的，这说明即使改进了奇异向量计算方案，也不会影响位温扰动结构。

图1 台风利奇马前 5个 TCSV（放大1000倍）对应的能量分布特征 （a1—a3、b1—b3.原方法和新方法对应的TCSV；a1、b1.不同TCSV对应的动能和位能分布，a2、b2.5个TCSV动能和位能及总能量均值的垂直分布，a3、b3.5个TCSV动能和位能及总能量的能量谱均值分布；KE：动能，PE：位能，TE：总能量；初始时刻为2019年8月8日12时）Fig.1 Energy distribution characteristics for the first five TCSVs （multiplied by 1000） corresponding to typhoon Lekima（a1—a3， b1—b3 respectively represent the TCSVs corresponding to the original SV computation method and the improved computation scheme of SV based on the hydrostatic equilibrium； a1，b1.the Kinetic energy （KE） and Potential energy （PE） distribution for different TCSVs，a2，b2.vertical distributions of the averaged KE，PE and Total energy （TE） for the five TCSVs，a3，b3.averaged energy spectra of KE，PE and TE for the five TCSVs； The initial time is 12:00 UTC 8 August 2019）

图2 基于原方法 （a1—a3、b1—b3） 和新方法 （c1—c3、d1—d3） 得到的台风利奇马第 3 个 TCSV（放大 1000 倍）的位温扰动结构 （a1—d1.初始扰动，a2—d2.TCSV 线性演变 48 h 后的扰动，a3—d3.TCSV 非线性演变 48 h 后的扰动；a1—a3 和 c1—c3.模式第26层的水平结构，b1—b3和d1—d3.沿25°N的垂直剖面；初始时刻为2019年8月8日12时）Fig.2 Structures of the potential temperature perturbation component of the third TCSV （multiplied by 1000） corresponding to the typhoon Lekima from the original SV computation method （a1—a3，b1—b3） and the improved computation scheme（c1—c3，d1—d3） of SV based on the hydrostatic equilibrium （a1—d1.TCSV at the initial time，a2—d2.TCSV after 48 h of linear evolving，a3—d3.TCSV after 48 h of nonlinear evolving；a1—a3，c1—c3.horizontal structures at the model level 26，b1—b3，d1—d3.vertical structures at 25°N；The initial time is 12:00 UTC 8 August 2019）

图3给出原方法和新方法台风利奇马对应的第3个TCSV（放大1000倍）的气压扰动分量在模式第26层的水平结构和在25°N的垂直结构。结果表明，初始时刻TCSV的气压扰动分量具有较大的差异，原方法对应的TCSV气压扰动分量的分布存在着强的局地性，表现为非常小尺度的水平分布结构，这种结构不利于TCSV非线性积分时的计算稳定性。奇异向量计算方法改进后，TCSV的气压扰动分量的大值区位于台风利奇马附近，更好地描述了热带气旋初值的不确定性，与位温扰动分量的分布更加协调，有利于TCSV非线性积分的稳定性。对于两种方法中气压扰动分量的线性演变和非线性演变特征，图3中呈现出类似的结构，这可能是因为与奇异向量的位温扰动、纬向风和经向风扰动量值相比，气压扰动分量的值过小（约小103倍），由于新方法对位温扰动、纬向风和经向风扰动的量值和结构影响较小，这两个方法中这些主导扰动分量的结构类似，从而在演变过程中气压扰动分量的特征由其他变量的演变特征主导，使得两个方法中气压扰动分量的演变特征类似。

图3 同图2，但为气压扰动分量的水平结构Fig.3 The same as Fig.2 but for the pressure perturbation component

此外，对比两种方法的气压扰动分量和位温扰动分量的分布，可以看出新方法中这两个扰动量分布特征更为协调，为产生更为协调的初值扰动奠定了基础，可能有利于集合预报技巧的提高。

3.2 初值扰动结构分析

基于两种奇异向量求解方法，分别针对北半球副热带（30°—80°N）、南半球副热带（80°—30°S）和低纬度热带气旋区域（以热带气旋为中心的10个经纬度区域）求解SV和演化SV，进而将南北半球副热带SV、演化SV和TCSV通过高斯线性组合产生初值扰动。有关初值扰动构造的细节参见李晓莉等（2019a）和霍振华等（2020）。

图4和5分别给出了原方法和新方法第5个初值扰动在模式第40层（大约200 hPa）的水平分布和沿25°N的垂直剖面。可以看出，奇异向量求解方法改进前、后初值扰动的位温扰动分量和纬向风扰动分量结构类似；奇异向量求解方法改进前，气压扰动分量的分布更加局地化，奇异向量求解方法改进后，气压扰动分量的分布可以更好地描述副热带中高纬度地区的斜压不稳定和热带气旋的初值不确定，在模式中、低层扰动值减小，与位温扰动分量的分布更加协调，可能有利于集合预报成员的积分稳定和集合预报技巧的提高。

图4 基于原 SV 求解方法得到的第 5 个初值扰动在模式第 40 层的水平分布 （a1—a3） 和沿 25°N 的垂直剖面 （b1—b3）（a1、b1.位温扰动，a2、b2.气压扰动，a3、b3.纬向风扰动；初始时刻为 2019 年 8 月 8 日 12 时）Fig.4 Horizontal structure at model level 40 （a1—a3） and vertical profile at 25°N （b1—b3） of the fifth initial perturbation based on the original SV computation method （a1，b1.potential perturbation，a2，b2.pressure perturbation，a3，b3.zonal wind perturbation；The initial time is 12:00 UTC 8 August 2019）

续图4Fig.4 Continued

图5 同图4，但SV求解方法为静力平衡奇异向量改进方法Fig.5 The same as Fig.4 but with the improved computation scheme of SVs based on the hydrostatic equilibrium

续图5Fig.5 Continued

4 静力平衡奇异向量改进方法对集合预报技巧的影响

基于两种奇异向量求解方法所获得的分目标区的SV，分别构建集合预报扰动初值，并开展了两组集合预报试验（2019年 8月 4、6、8、11、13和15日每天进行，12时起报），分析新方法对集合预报技巧的影响。基于所有个例的集合预报技巧评分均值分析了新方法对北半球和南半球等压面要素集合预报、热带气旋路径集合预报和中国地区24 h累计降水概率预报技巧的影响。由于新方法对北半球和南半球等压面要素集合预报的影响类似，在此仅给出北半球等压面要素集合预报技巧结果，南半球结果不再赘述。

4.1 北半球等压面要素集合预报技巧

4.1.1 均方根误差和集合离散度

均方根误差（Root Mean Square Error，RMSE）度量的是集合平均值与真值（观测值或分析场，本试验使用的是业务分析场作为真值）的偏差程度，均方根误差越小，则预报越接近真值。集合离散度（SPREAD）衡量的是集合成员与集合平均值的差异，离散度过小，意味着不够发散，不足以涵盖到真值，容易导致漏报；离散度过高，则意味着集合成员太过发散，容易导致虚警。

图6给出了北半球等压面要素均方根误差和集合离散度随预报时效的演变。结果表明，奇异向量求解方法改进后，预报前期，集合平均RMSE变化不大，集合离散度增大，集合离散度和集合平均RMSE的关系改善；预报后期，集合平均RMSE和集合离散度变化不大。

4.1.2 Outlier评分

预报失误概率Outlier评分是由Talagrand直方图进一步发展出的评分，是分析（观测）值落在Talagrand直方图中最两端的两个盒子（bins）中的频率，即分析（观测）落在集合成员的预报值区间之外的平均频率，以表示集合预报系统的漏报率。一个具有n个成员的理想集合预报系统，其实况落在集合预报区间之外的概率期望值为2/（n+1）。Outlier评分是评估集合预报系统可靠性指标之一，Outlier值偏离2/（n+1）越多，说明预报失误的概率越大；Outlier值偏离 2/（n+1）越少，则漏报率越低，系统越可靠。本研究中，集合成员个数为31，从而Outlier值的理想值是6.24%。

图7分别给出了北半球等压面要素Outlier评分随预报时效的演变。结果表明，奇异向量计算方法改进后，从整个预报时段来看，GRAPES全球集合预报系统等压面要素的Outlier评分与概率期望值更接近，即预报失误概率减小，概率预报技巧增加。

4.1.3 CRPS 评分

连续分级概率评分（Continuous Ranked Probability Score，CRPS）是评估集合预报系统整体性能常用的评分指标，代表观测和预报的累积分布函数的差别。CRPS值越大，观测和预报的累积密度函数的差别越大，表示集合预报系统的预报能力越低。

图8分别给出了北半球等压面要素CRPS评分随预报时效的演变。结果表明，奇异向量计算方法改进后，GRAPES全球集合预报系统的北半球等压面要素CRPS评分有一定降低，即CRPS概率预报技巧有所提高。

图6 北半球等压面要素的均方根误差和集合离散度随预报时效的演变 （黑色实线：对照预报的均方根误差，蓝色实线 （虚线）：基于原奇异向量求解方法得到的集合平均预报的均方根误差 （集合离散度），红色实线 （虚线）：用静力平衡奇异向量改进方法得到集合平均预报的均方根误差（集合离散度）；a.500 hPa 位势高度，b.850 hPa 温度，c.250 hPa 纬向风速，d.250 hPa 经向风速）Fig.6 Evolutions of RMSE and SPREAD for the northern hemisphere isobaric variables （black solid line: RMSE of the control forecast，bule solid （dashed） line: RMSE （SPREAD） of the ensemble mean （EM） obtained with original SVs computation method，red solid（dashed） line: RMSE （SPREAD） of EM obtained with the improved computation scheme of SVs based on hydrostatic equilibrium；a.geopential height at 500 hPa，b.temperature at 850 hPa，c.zonal wind at 250 hPa，d.meridional wind at 250 hPa）

图7 北半球等压面要素Outlier评分随预报时效的演变 （蓝色实线：原SV求解方法对应的结果，红色实线:静力平衡奇异向量改进方法对应的结果，黑色虚线：Outlier期望值；a.500 hPa 位势高度，b.850 hPa 温度，c.250 hPa 纬向风，d.250 hPa 经向风）Fig.7 Evolution of the Outlier scores for the northern hemisphere isobaric variables （blue solid line: the original SV computation method，red solid line: the improved computation scheme of SV based on hydrostatic equilibrium，black dashed line: the Outlier expectation；a.geopotentail height at 500 hPa，b.temperature at 850 hPa，c.zonal wind at 250 hPa，d.meridional wind at 250 hPa）

续图7Fig.7 Continued

图8 北半球等压面要素CRPS评分随预报时效的演变 （蓝色实线：奇异向量求解方法结果，红色实线：静力平衡奇异向量改进方法结果；a.500 hPa 位势高度，b.850 hPa 温度，c.250 hPa 纬向风，d.250 hPa 经向风）Fig.8 Evolutions of the CRPS scores for the Northern isobaric variables （blue solid line: the original SV computation method，red solid line: the improved computation scheme of SV based on the hydrostatic equilibrium；a.geopotentail height at 500 hPa，b.temperature at 850 hPa，c.zonal wind at 250 hPa，d.meridional wind at 250 hPa）

4.2 台风集合预报技巧

图9a给出了所有台风个例路径预报误差和集合离散度均值随时间的演变。图9b和c分别给出了路径误差和离散度的箱线图。从图9a来看，相对于原方法试验，新方法试验中台风路径集合离散度增大，集合平均路径预报误差在预报前期略微增大，预报后期有一定减小，集合平均路径预报误差和路径集合离散度的关系得到改善。从箱线图（图9b、c）来看，不同台风集合预报技巧随时间变化各有不同，随着预报时间延长，对照预报、集合平均误差和离散度呈现增大的趋势。从第一四分位数和第三四分位数来看，6—60 h集合平均误差和对照预报误差相差不大，66—120 h新方法集合平均误差小于对照预报和旧方法集合平均。在整个预报时段，无论从第一四分位数、中位数还是第三四分位数来看，奇异向量计算方法改进之后，台风路径离散度明显增大，与集合平均路径预报误差更接近。对原方法和新方法试验中的集合平均路径误差（集合离散度）进行了配对t检验，结果表明，对6—120 h内的任一预报时次，原方法和新方法试验集合平均路径误差的差异均没有通过95%显著性检验，原方法和新方法试验集合离散度的差异均通过了95%的显著性检验。这说明，新方法显著地增加了台风路径集合离散度，对集合平均路径预报误差没有显著影响。

图9 台风路径预报误差和路径集合离散度随预报时效的演变 （Control：台风路径预报误差，EM_Old （EM_New）：原 （新） SV求解方法台风集合平均路径预报误差，SPD_Old（SPD_New）：原 （新） SV求解方法台风路径集合离散度；a.所有台风个例对应的路径预报误差和集合离散度均值，b.台风路径误差箱线图，c.台风路径集合离散度箱线图）Fig.9 Evolutions of the yphoon track forecast error and ensemble spread （control：the typhoon track forecast error of control forecast，EM_Old （EM_New）：the ensemble mean track forecast error corresponding to the original （new） SVs computation method，SPD_Old （SPD_New）：the typhoon track ensemble spread corresponding to the original （new） SVs computation method；a.averaged typhoon track forecast error and ensemble spread for all TC cases，b.boxplot for the track forecast error，c.boxplot for the track ensemble spread）

4.3 中国地区24 h累计降水概率预报技巧

AROC评分和BS评分是业务上检验集合预报性能常用的评分，用来检验分级降水的预报能力。AROC评分越大，BS评分越小，则降水概率预报技巧越高。图10和11分别给出了中国地区小雨、中雨、大雨、暴雨各量级24 h累计降水概率预报技巧评分AROC和BS随预报时效的演变。结果表明，奇异向量计算方法改进后，对于小雨、中雨、大雨，AROC评分增大，BS评分减小，即中国地区24 h累计降水概率预报技巧提高；对于暴雨，AROC评分某些时次增大，某些时次减小，BS评分在整个预报时段内均减小。总体来看，采用新方法的集合预报能提高中国地区暴雨24 h累计降水概率预报技巧。

图10 中国地区小雨、中雨、大雨、暴雨各量级24 h累计降水概率预报技巧评分AROC随预报时效的演变（a.小雨，b.中雨，c.大雨，d.暴雨）Fig.10 Evolution of the precipitation probability forecast skill score AROC for 24 h accumulated precipitation over China（a.light rain，b.moderate rain，c.heavy rain，d.rainstorm）

5 总结和讨论

基于奇异向量的初值扰动方法已应用于GRAPES全球集合预报系统，为解决部分奇异向量非线性积分崩溃的问题，基于原有的GRAEPS全球奇异向量计算方案，发展了静力平衡奇异向量改进方法，在求解奇异向量时，将初始时刻的气压扰动分量根据静力平衡关系得到，改进了原有的奇异向量求解方法。基于两种奇异向量求解方法构建集合预报扰动初值，开展了集合预报试验，分析了静力平衡奇异向量改进方法对集合预报技巧的影响。研究结果表明：

图11 同图10，但为 BS 评分Fig.11 The same as Fig.10 but for Brier score

（1）静力平衡奇异向量改进方法通过产生协调的气压扰动和位温扰动场，解决了奇异向量非线性积分崩溃的问题，消除了原来不利于积分稳定性的气压扰动过于局地化的小尺度结构，使得气压扰动分量的分布更广，与位温扰动分量的分布更协调。同时，静力平衡奇异向量改进方法对奇异向量中位温扰动分量和纬向风扰动分量结构影响较小，显示了静力平衡奇异向量改进方法的合理性。

（2）基于新奇异向量求解方法分目标区的SV构建出的集合预报初值扰动中气压扰动结构与其直接计算出奇异向量中的气压扰动分量的特征类似，明显减小了气压扰动值的小局地化扰动特征，与初值扰动中位温扰动分量的分布更加协调。

（3）从集合预报试验的结果来看，新奇异向量求解方法在预报前期对集合平均的均方根误差影响不大，但一定程度增加了集合离散度，使得集合离散度和集合平均的均方根误差的关系得到改善；预报后期，集合平均的均方根误差和集合离散度变化不大。从整个预报时段来看，基于新奇异向量求解方法的集合预报试验中Outlier评分与概率期望值更接近，即预报失误概率减小，同时，CRPS概率预报技巧有所提高。

（4）基于新奇异向量求解方法能改进GRAPES全球集合预报的台风路径集合离散度，使得集合平均路径预报误差和路径集合离散度的关系得到改善，并提高了中国地区小雨、中雨、大雨、暴雨各量级24 h累计降水概率预报技巧。

总之，新发展的静力平衡奇异向量改进方法不仅解决了奇异向量非线性积分崩溃的问题，也对GRAPES全球集合预报效果有正贡献，下一步，会将此方法应用于GRAPES全球集合预报业务系统。需要指出的是，文中初步采用静力平衡奇异向量改进方法解决了采用湿线性物理方案时TCSV计算的稳定性问题，随着GRAPES切线性模式中线性化物理过程的完善，在未来的研究中，将继续深入开展GRAPES全球模式湿线性物理过程对TCSV及热带气旋集合预报技巧的影响研究。此外，尽管奇异向量方法已成功应用于集合预报，但奇异向量基于线性理论，不能反映非线性物理过程对天气和气候可预报性的影响，在研究非线性模式中有限振幅初始扰动引起的可预报性问题方面具有局限性（Gilmour，et al，1997；Anderson，1997）。正交条件非 线性最 优 扰动 （Conditional Nonlinear Optimal Perturbations，CNOP；Duan，et al，2016）方法可以获得非线性模式中增长较快的初值扰动，在一定程度上克服了奇异向量方法的线性局限性。CNOP在集合预报中的应用研究（Duan，et al，2016；Huo，et al，2019a，2019b）表明，使用正交 CNOP 方法可以得到较高的集合预报技巧，可以有效提高台风路径的预报技巧，是一种有潜力的集合预报方法。为弥补奇异向量方法在GRAPES全球集合预报初值扰动构造方面的线性局限性，在未来的研究中，将开展正交CNOP方法在GRAPES全球集合预报中的应用研究，以期提高GRAPES全球集合预报尤其是台风集合预报方面的技巧。