1 引　言

天气模式为了能够使用较长的时间步长，通常使用半隐式时间积分算法，对低频波动，如缓慢变化的天气尺度波动显式积分；对高频波动，如快速变化的重力波、声波等隐式求解。在半隐式时间积分算法的求解过程中，一般需要引入参考态，将温度、气压等预报变量分解为参考态和扰动之和，并对预报方程进行线性化分离，将其分解为线性项和非线性项，参考态的相关计算一般划入线性项，扰动划入非线性项。为了保证半隐式时间积分的精度，该类方法要求线性项明显大于非线性项。因此，参考态的选取非常重要，如果参考大气可以尽量接近模式大气，就可以减小扰动量，提高半拉格朗日空间插值的计算精度，同时减小非线性项，提高半隐式时间积分的精度。另外，陡峭地形处的虚假响应和平流计算一直是限制半拉格朗日模式时间步长选取的主要难题，对积分稳定性也有显著影响。Ritchie等（1996）研究指出，地形的强迫作用在积分过程中不随时间发生变化，可以通过参考态将其扣除，以此来提高陡峭地形处的气压梯度力和半拉格朗日平流过程的计算精度。

前期的研究工作（苏勇等，2018）在调研不同参考大气的设计方法之后，在GRAPES_GFS模式（陈德辉等，2006；沈学顺等，2011，2017；王金成等，2017；张华等，2004；刘艳等，2019）中引入了不随时间变化且满足静力平衡的三维参考大气，替换由等温大气构造的一维参考大气，使参考大气在积分过程中可以尽量接近模式大气。苏勇等（2018）重新推导了引入三维参考大气之后动力框架的求解过程，并通过一系列的理想试验测试了理论方法以及代码实现的正确性，验证了三维参考大气在GRAPES_GFS中可以有效地提高框架的计算精度。但前期的若干个理想试验中都是选取试验初值为其三维参考态，初始扰动为0，该设置仅适用于理想试验。实际资料预报的情况下不能直接将初值设置为三维参考态，需要根据动力框架的求解方法，考虑精度、稳定性等各个方面的因素，构造合理的、适用的三维参考态。

本研究工作主要解决以下两个问题：（1）实际资料预报情况下三维参考大气的设计，也就是如何根据实际大气给定其三维参考态；（2）评估GRAPES_GFS引入三维参考大气之后实际资料预报效果的改进，对业务系统预报能力的影响。

早期的研究工作主要集中在利用等温大气构造的参考态，或者温度仅随高度变化的一维参考廓线。曾庆存（1963）在ECMWF（European Centre for Medium range Weather Forecast）谱模式的静力学方程组中通过扣除温度和高度的参考态来提高预报精度。Chen等（1990）在ECMWF的业务模式中通过静力扣除参考大气的方法提升了模式的整体预报性能。Zhang等（1990）将Chen等（1990）的方法进行了扩展，在RMRC（Bureau of Meteorology Research Center）的全球谱模式中引入了随高度和纬度变化的二维参考大气，在天气学模拟中得到了正的预报效果。Wu等（2008）在NCAR（National Center for Atmospheric Research）的CAM3（Community Atmosphere Model version 3）中，通过温度方程中扣除随高度变化的一维廓线，地面气压预报方程中扣除和地形相关的地面气压分量，提高了区域和全球气候尺度模拟的效果。

近期的研究工作中，ECMWF的IFS（Integrate Forecast System）和英国气象局ENDGame（Even Newer Dynamics for General atmospheric modelling of the environment）最具有先进性和代表性。IFS模式质量方程和能量方程的预报变量使用地面气压和温度，Temperton等（2001）从其质量方程中扣除和地形相关的地面气压分量，然后构造一个和该地面气压分量以及模式垂直层次共同相关的温度三维参考态，并从能量方程中将其扣除。该方法有效地提高了IFS模式中陡峭地形处半拉格朗日平流和气压梯度力的计算精度。英国气象局的做法（Wood，et al，2014）是目前理论上最为先进的方案，将上一时步的三维场作为参考态，每一步的扰动态是相对于上一步的变化量，相当于采用随时间变化的三维参考大气，可以使参考大气在积分过程中始终非常接近模式大气。在具体计算过程中，方程中需要包含参考态的时间变化项，每个积分步开始时需要在每个点的垂直方向上调整位温的参考态使其满足单调递增，然后根据静力平衡关系导出气压的参考态。对于参考态的水平偏导数项，英国气象局采用方程两端同时增减自己设计的垂直运动项，巧妙地降低了计算的复杂性。但整体而言，英国气象局关于三维参考大气的求解过程较为复杂，且计算代价较大，其他模式中不易实现。

对比ECMWF和英国气象局的方法，可以看到两者侧重点不同：ECMWF主要关注陡峭地形处预报效果的提升，根据文献中的公式进行计算，其在地面气压和温度预报方程中扣除地形相关的参考态部分，约占全量的40%，参考态和地形高度相关，但并不贴近实际大气；英国气象局的方法使参考态尽可能的贴近实际大气，根据GRAPES_GFS每时步温度的变化进行估计，其参考态约占全量的99%以上，扰动态不足1%。虽然英国气象局没有直接根据地形构造参考态，但参考态中实际上包含了各个尺度的天气系统以及地形相关的信息。英国气象局的方法大幅度减小扰动量的数量级，可以使非线性项远小于线性项，有助于进一步提高半隐式时间积分算法的精度。

2 参考态构造

如何根据模式大气给定三维参考大气的参考态，是本研究的重点。初期尝试在GRAPES_GFS中参考IFS的处理方法，先用地面气压计算出地面的无量纲气压，然后根据等温大气或者随高度变化的温度廓线，利用静力平衡关系逐层积分得到三维参考大气。但将该三维参考态放入模式中进行实际预报试验时，动力框架在求解亥姆霍兹方程时无法收敛，模式很快积分溢出。如果对参考大气进行空间平滑，则积分初始若干步可以收敛，其后仍无法收敛最终溢出。尝试对问题进行各个方面的分析，但由于GRAPES_GFS和IFS两个模式的方程组设计、框架求解过程等存在很大差异，地形扣除的方法无法一致，暂未查明出问题的原因。

广州热带气象研究所在其高分辨率版本的GRAPES有限区域模式中使用和文中相似的三维参考大气，在实际资料预报中采用对初值进行9点平滑作为参考大气的方法。本部分的研究工作借鉴其思路，在GRAPES_GFS系统循环预报的情况下，采用如下流程构造三维参考大气和模式初值：（1）从同化循环的初值中读取无量纲气压（ $\varPi $ ）；（2）对 $\varPi $ 进行9点平滑得到无量纲气压参考态的预估值（ ${\overline {\varPi} ^ * }$ ）；（3）根据静力平衡关系，由 ${\overline {\varPi} ^ * }$ 计算位温参考态的预估值（ ${\overline \theta ^ * }$ ）；（4）在每个格点的垂直方向上利用冒泡法对 ${\overline \theta ^ * }$ 进行排序，得到位温参考态（ $\overline \theta $ ），确保 $\overline \theta $ 在垂直方向单调递增；（5）根据静力平衡关系，由 $\overline \theta $ 重新计算得到无量纲气压参考态（ $\overline {\varPi} $ ）；（6）计算无量纲气压扰动量： $\varPi ' = \varPi - \overline \varPi $ ；（7）根据静力平衡关系，由 $\varPi '$ 计算位温扰动量（ $\theta '$ ）；（8）计算位温： $\theta = \overline \theta + \theta '$ 。

本研究将该方法称为初值光滑法，该方法的优点是：（1）初始场中参考大气非常接近模式大气；（2） $\overline \theta $ 在垂直方向单调递增确保参考大气的静力稳定，利于模式稳定积分；（3）初值中全量的 $\varPi $ 和 $\theta$ 也保持严格的静力平衡关系，初值中没有垂直速度，利于模式的稳定积分。但该方法也有缺陷，方法中引入9点平滑的目的是使参考大气较为光滑，以保证陡峭地形处参考态偏导数的光滑性，利于模式的稳定积分。但后期的工作中显示，对参考态直接进行空间平滑会导致陡峭地形处参考态和地形在空间分布上不一致，陡峭地形处模式会计算出与地形分布不匹配的垂直速度，进而导致不合理的、类似点阵状的降水分布。

借鉴以上初值光滑方法的优点和缺点，本研究设计了气候平均的方法，仍以GRAPES_GFS系统循环预报情况为例，具体流程如下：（1）用同化循环的初值求月平均的无量纲气压，作为 ${\overline {\varPi} ^ * }$ ；（2）后续流程同初值平滑方法的第（3）—（8）步。

该气候平均方法避免了初值光滑法中的空间平滑导致的陡峭地形处参考态和地形分布不一致的问题，通过时间平滑的方法构造三维参考大气，保留了初值光滑方法的优点（2）和（3），相对来说，扰动量略大于初值光滑方法，但仍明显小于一维参考大气。该方法中月平均态的三维参考大气较为光滑，可以保证参考态的水平和垂直偏导数等计算中间变量在高纬度或者复杂地形处的光滑度，积分非常稳定。另外，该方法对时间平均的跨度并不敏感，如采用月平均场、季节平均场可以得到基本一致的预报结果，因为全量并未发生变化，只是很少的一部分信息在参考态和扰动之间相互转移。

本研究中GRAPES_GFS系统都采用中国气象局数值预报中心发布的2.4版本，采用该版本进行循环同化，任选循环过程中一个时次做预报。对比三维参考大气和一维参考大气在积分1 d之后扰动量沿180°经线的分布（图1），可以看出，一维参考大气扰动量的量级很大， $\varPi '$ 的数值范围大约为0.01—0.08， $\theta '$ 大约为40—100 K，对流层高层达到100 K，约占全量的25%。采用初值平滑和气候平均两种方法构造三维参考大气之后，扰动量相对于一维参考大气缩小了一个量级， $\varPi '$ 为−0.002—0.004， $\theta '$ 为−15—15 K。对比初值平滑和气候平均两种方法，扰动量的量级和分布形态较为接近，无明显区别。

进一步对比初值平滑和气候平均两种方法，图2给出模式积分1 d之后 $\overline \theta $ 沿东西向水平偏导数在南半球的分布。因为在循环预报的过程中，模式的初值在高纬度地区带有一定的计算噪音，平滑或者滤波算法无法彻底消除这些噪音，再加上经、纬度网格在高纬度区域水平格距较小，所以在对 $\overline \theta $ 求水平偏导数之后，计算噪音更加明显。初值光滑方法的情况（图2a）下， $\overline \theta $ 的水平偏导数在南极大陆，特别是陡峭地形附近，正、负相间分布的计算模态非常明显。参考态的水平、垂直偏导数等在模式动力框架求解中非常重要，且不会随积分时间得到修正，上述计算噪音会影响高纬度区域的计算精度，在某些情况下也会造成积分溢出。气候平均方法的情况（图2b）下， $\overline \theta $ 的水平偏导数相对来说比较光滑，只有在陡峭地形附近梯度较大，没有明显的计算模态。对比个例试验的预报效果也可以看出，气候平均方法构造的三维参考大气，在高纬度地区预报效果较好，且在一年的冷启动预报试验和两个月的循环同化预报试验中都没有积分溢出现象，所以本研究最终选取气候平均方法构造三维参考大气。

对比本研究中的气候平均法和英国气象局三维参考态构造方法可以看到，两者存在相似之处，后续可以逐步地将气候平均法向英国气象局的方法进行调整。气候平均法中采用月平均场作为三维参考态，参考态满足静力平衡、位温单调递增且不随积分时间发生变化；英国气象局的方法用上一时步的场作为三维参考态，参考态同样满足静力平衡、位温单调递增，但参考态随着模式积分每步都发生变化。后续工作中考虑在GRAPES_GFS中引入参考态随时间的变化项，并将气候平均场改为模式初值，这样气候平均法就可以调整为英国气象局的方法，但由于两个模式方程组设计以及求解方法的差异，如何保证模式积分稳定性会是该部分工作的难点。

3 实际资料预报试验 3.1 动能谱

将数值模式的动能谱与实际大气动能谱进行比对分析，是评估检验模式框架能量转换、摩擦耗散等性能的有效途径（郑永俊等，2008）。实际大气动能在各个波段随波数衰减的统计规律，在大尺度范围服从 ${k^{ - 3}}$ 的斜率，在中尺度范围服从 ${k^{ - 5/3}}$ 的斜率。利用标准化的动能谱的计算程序，对GRAPES_GFS采用一维参考大气和三维参考大气两种情况下的结果进行诊断，并与NCEP（National Centers for Environmental Prediction）准业务系统FV3gfs（Finite-Volume 3 GFS system）模式的结果进行对比。其中，GRAPES_GFS模式水平分辨率0.25°，垂直60层，时间步长300 s，初始场采用NCEP的FNL（Final operational global analysis data）资料，冷启动预报8 d，预报结果为标准等压面上的经纬度网格数据；FV3gfs模式的预报结果来自于国家气象信息中心下载的FV3gfs平行试验数据，为标准等压面上的经纬度网格数据，水平分辨率0.5°。FV3gfs模式原始格点为立方球网格，水平分辨率约13 km，垂直62层。因为两个模式的数据在分辨率等方面存在差异，所以FV3gfs模式的动能谱只用来做参考性的对比。3组试验均为5个样本，试验日期为2018年7月1日、8月1日、9月1日、10月1日、11月1日，起报时间为12时（世界时，下同），对积分24和72 h的动能谱分别进行对比，先计算每个样本的动能谱，然后平均得到5个样本的平均动能谱。

对比过程中发现，不同个例、不同预报时长之间动能谱的分布差异都不大，所以此处仅选择积分24 h的平均动能谱进行对比。GRAPES_GFS和FV3gfs在长波范围内动能随波数的衰减规律基本一致，在中短波范围内，二者的衰减率相对于参考斜率都有一定偏离，但相对来说，FV3gfs更接近参考线，在中小尺度范围内，能量衰减较小，与实际大气更为相符。对比GRAPES_GFS中采用一维参考大气和三维参考大气的两组试验，在100 hPa上三维参考大气的能量衰减明显小于一维参考大气，在中小尺度范围内可以拟合出 ${k^{ - 5/3}}$ 的关系，动能谱的分布更接近实际大气（图3）。

3.2 四维变分循环试验预报效果

基于GRAPES_GFS 2.4版本，加入三维参考大气相关程序，其中三维参考大气的构造方法选择气候平均法。在资料同化科相关同事的协助下，进一步修改了四维变分（Zhang，et al，2019；刘永柱，2017）模块中切线性部分和伴随部分的代码，之后进行了四维变分循环试验，并与采用一维参考大气的控制试验进行了对比。四维变分循环试验的时段为2016年6月1日到7月31日，对每天12时的分析场做8 d预报，水平分辨率0.25°，垂直60层，两组试验的时间步长一致，均为300 s。两组试验在同化循环和模式预报的过程中都可以保证稳定积分，无溢出现象。采用三维参考大气增加了框架求解过程的复杂性，但对模式计算效率的影响较小，模式积分耗时增加了2%—3%。

图4为四维变分循环试验预报试验的综合评分卡，分东亚、北半球、南半球、赤道4个区域分别对风场、温度场、高度场进行检验，红色表示正效果，绿色表示负效果，灰色为中性，箭头越大效果越显著。可以看到采用三维参考大气之后，对于夏季两个月来说，模式的综合预报性能得到了提升，北半球和东亚区各个预报变量前3 d或者前5 d范围内预报效果都有所改善，南半球预报效果也略有提升，赤道地区风场有所改善，温度、高度场在个别层次上预报性能下降。

进一步对检验结果进行分析，从北半球和南半球500 hPa高度场的距平相关系数（图5）可以看出，使用三维参考大气之后高度场的距平相关系数并没有得到明显的提升，说明槽脊的分布形态并没有发生显著的变化，但从高度场的预报偏差（图6）中可以看到，三维参考大气明显地减小了一维参考大气在50 hPa以下的负预报偏差，对流层的高度场模拟更为合理。对于温度场的检验也有类似的结果。从风场的预报偏差（图7）中可以看到，采用三维参考大气之后，显著地减小了对流层顶部100 hPa附近的预报偏差，对流层中下部整体偏差结构变化不大，偏差数值略有减小。

GRAPES_GFS在长期积分过程中存在质量损失较为严重的问题（苏勇等，2016），采用一维参考大气的情况（图8实线）下，对于两个月的批量试验来说，全球平均的海平面气压下降了约0.7 hPa，相比于国际上先进预报中心的业务模式来说（Diamantakis，2013），长期积分质量守恒性能较差。在采用三维参考大气的情况（图8虚线）下，全球平均的海平面气压下降了约0.2 hPa，大幅度提升了模式在长期积分过程中的质量守恒性。

4 结论和讨论

（1）对于实际资料预报情况下的三维参考态的给定方法，对比分析了ECMWF的IFS中地形扣除的方法以及英国气象局ENDGame中参考态随时间变化的方法，并借鉴广州热带气象研究所初值光滑方法的优缺点，结合GRAPES_GFS的现状，采用气候平均方法，用月平均场作为三维参考大气，使扰动量的数量级相对于一维参考大气明显减小，同时可以保证参考态较为光滑，确保模式的稳定积分。

（2）通过对模式动能谱的比对分析，可以看到使用三维参考大气之后，模式在高层的能量耗散明显减小，可以拟合出更为合理的能量衰减规律。

（3）通过夏季两个月的循环预报试验，可以看到使用三维参考大气之后，模式的综合预报性能有所提升，对流层高度场、温度场的预报偏差有所减小，长期积分过程中模式质量损失较为严重的问题也得到了明显的缓解。

后续工作中，以下几点问题值得讨论：

（1）构造参考大气的基本考虑是使参考态尽量接近模式大气，尽可能的减小扰动量的数量级，如英国气象局的做法（Wood，et al，2014），每一步求解过程中的参考态是上一步的变量场，相当于随时间变化的三维参考大气。在GRAPES_GFS中，基于本研究的三维参考大气，如果加入参考态的时间变化项，也可以实现随时间变化的三维参考大气。但在这种情况下，参考态中各个尺度的信息丰富，和地形的相关度也很高，如何保证参考态的计算精度，确保模式积分的稳定性，是需要重点解决的问题。

（2）GRAPES_GFS现在垂直方向60层，模式层顶高度为36 km，约4 hPa，近期工作正在将垂直层次增加至87层，模式层顶高度设置为63 km，约0.1 hPa。在抬高模式层顶的情况下，对于基于等温大气构造出的一维参考大气来说，在平流层的上部和中层，位温的扰动态将和参考态量级相当，线性化分离的求解方案变得没有意义，采用三维参考大气则可以解决这一问题，在平流层和中层也可以保证扰动态明显地小于参考态，确保时间、空间积分算法的精度。

（3）气候模式的动力框架在求解过程中通常也需要引入参考大气，将预报变量分解为参考态和扰动态。Wu等（2008）通过改进CAM3中温度随高度变化的一维参考廓线，提高了区域和全球气候模拟的预报效果。本研究中的三维参考大气方法，在气候模式中也有一定的应用前景，如在气候模式的积分过程中，将年平均或者季节平均的气候态作为参考态扣除，就可以大幅度的减小扰动态的数值，同时扣除地形强迫作用对平流计算的部分影响，预计可以有效地提高气候模式动力框架的空间、时间离散化精度。