不同观测误差确定方法对地基GNSS-ZTD 资料同化预报效果影...

奥鹏网院作业

不同观测误差确定方法对地基GNSS-ZTD 资料同化预报效果影响的对比分析*
公 颖1 李得勤1 仲跻芹2 周晓珊1 杨 森1 潘 晓1

1.中国气象局沈阳大气环境研究所，沈阳，110166

2.中国气象局北京城市气象研究院，北京，100089

摘 要 利用2016 年6—8 月华北—东北地区的地基全球卫星导航系统的天顶总延迟（GNSS-ZTD）观测资料、东北区域中尺度数值预报系统，以2016 年6—8 月的13 d 强降水为例，开展基于Desroziers 等（2005）理论的Des 方法和传统方法进行观测误差确定的天顶总延迟资料同化对比试验研究，探讨Des 方法相对于传统观测误差确定方法对天顶总延迟资料同化预报效果的影响，并以未做天顶总延迟资料同化的试验为对照试验，考察天顶总延迟资料在数值模式中的同化应用效果。结果表明：（1）Des 方法得到的天顶总延迟观测误差诊断值较为合理，诊断值站点间差别较大，说明逐站进行观测误差诊断的必要性；（2）天顶总延迟资料同化使强降水的强度、落区预报性能得到提高，使温、湿、风等要素的预报与观测接近，Des 方案同化分析、预报效果优于传统方案；（3）对2016 年7 月25 日华北—东北强降水过程进行了同化预报分析，整体而言，天顶总延迟资料同化有效增强了对流层中低层初始湿度场，修正了积分初期水凝物含量与位置，进而改善了降水预报效果，修正了对照试验对辽宁东部地区强降水的明显漏报，且通过降水的反馈作用改进了温度与风场预报效果。基于Des 方法逐站诊断观测误差相比传统方法得到的观测误差更为合理，因此能够提高天顶总延迟资料的同化预报效果，同化天顶总延迟资料能够提高降水及温、湿、风等气象要素的预报水平。

关键词 观测误差确定，地基全球卫星导航系统的天顶总延迟（GNSS-ZTD）数据，资料同化，三维变分

1 引言
使模式初值能够更加准确地描述大气湿度分布进而改善降水预报水平一直是资料同化技术致力于解决的重要问题之一。地基全球卫星导航系统（GNSS）反演的大气水汽信息可被同化进数值模式，以提高初值场湿度要素分布的准确性（Vedel，et al，2001；Poli，et al，2007；陈敏等，2010；仲跻芹等，2017；Lindskog，et al，2017；Mascitelli，et al，2019；Rohm，et al，2019）。可同化进数值预报模式的GNSS 数据有天顶总延迟（ZTD）和大气可降水量（PW）两种，由ZTD 计算大气可降水量需要假定各向延迟同性，需要气压和温度等气象观测资料，这些都增加了大气可降水量的计算误差（丁金才，2009），而直接同化天顶总延迟的误差来源比同化大气可降水量少，因此，地基GNSS 观测资料在数值预报中的同化应用逐渐更多地集中于天顶总延迟资料（Vedel，et al，2001；Poli，et al，2007；Bennitt，et al，2012；Arriola，et al，2016；Yang，et al，2020；Giannaros，et al，2020）。

观测误差的确定是目前天顶总延迟资料同化工作中的关键与难点之一（Poli，et al，2007；Macpherson，et al，2008 ；Boniface，et al，2009；Bennitt，et al，2012；仲跻芹等，2017）。变分同化过程中，观测误差用于确定分析中观测资料的权重（Macpherson，et al，2008；Waller，et al，2016），是观测信息能否被充分有效利用以及分析场能否对大气状态进行最优估计的重要决定因素，由于观测真值无法获得，观测误差的确定成为当前资料同化中一个主要的困难与挑战（Chapnik，et al，2004）。到目前为止，绝大多数的天顶总延迟资料同化中，观测误差的确定方法均较为粗略（以下称传统方法）：通常基于天顶总延迟观测背景差（以下称Omb）标准差确定观测误差（Boniface，et al，2009；Yan，et al，2009a；Bennitt，et al，2012）；也有根据天顶总延迟平均值与Omb 标准差间的稳定拟合关系确定观测误差（Macpherson，et al，2008）；更多的是基于经验直接给定观测误差（Vedel，et al，2004；Poli，et al，2007；仲跻芹等，2017）。针对资料同化中普遍存在的观测误差不够精确的问题，有研究（Chapnik，et al，2006；王金成等，2015）表明，对观测误差进行诊断和调谐能够减小预报误差。

当前国际上有4 种具代表性的观测误差诊断理论与方法：（1）H-L 方法（Hollingsworth，et al，1986），该方法在观测误差空间不相关、预报误差空间相关等假设基础上，由新息增量样本计算和推算得到水平距离为0 的背景场误差协方差，在对应的新息增量协方差中将其扣除，即可获得观测误差协方差；此方法局限在探空观测站稠密的大陆地区使用（Kuo，et al，2004；王金成等，2015）。（2）最大似然方法（Dee，et al，1995，1999），该方法基于新息增量，通过最大似然估计确定误差协方差模式参数，进而得到观测误差协方差。（3）Desroziers 等（2001）基于观测目标函数最小值的期望值与实际值来计算观测误差调整参数；此方法仅适用于变分分析框架。（4）Desroziers 等（2005）提出基于Omb 与观测分析差（以下简称Oma）诊断观测误差协方差的方法，详见3.3.2 节；其优点是适用于任何分析框架。方法（4）是目前国际上较为流行的观测误差诊断方法，该方法应用于变分同化系统时需采用一次诊断的方式，而非理论上应采用的不动点迭代方法（Desroziers，et al，2005），因为在相关的观测误差被假定为不相关的变分同化系统中迭代方法不可行（Waller，et al，2016；Todling，2015），而这种一次诊断方法在业务天顶总延迟的三维变分（3DVar）同化中的适用性及其与传统观测误差确定方法的应用效果对比，尚未有过详细而系统的研究。

本研究利用2016 年6—8 月华北及辽宁区域的天顶总延迟观测资料，以东北区域中尺度数值预报系统为平台，开展了基于Desroziers 等（2005）理论的Des 方法和（普遍使用的）传统方法进行观测误差确定的天顶总延迟资料同化预报效果批量对比试验与检验，分析基于Des 方法逐站诊断观测误差相比传统方法对天顶总延迟资料同化预报效果的影响以及该资料的同化应用效果，以期能够为天顶总延迟资料同化以及资料同化中观测误差的诊断提供一些有益的参考。

2 东北区域中尺度数值预报系统简介
中国气象局沈阳大气环境研究所基于美国国家大气研究中心发展的WRF 模式和WRFDA 资料同化系统发展了东北区域中尺度数值预报系统，为东北区域的气象预报与服务提供了重要的技术支撑。该系统主模式为WRFv3.8.1，同化模块为WRFDAv3.8.1，采用3DVar 技术进行各种观测资料的同化，设计两层嵌套，空间分辨率分别约为9、3 km，分别覆盖全中国地区和华北东部—东北地区。所用物理过程和参数化方案包括：Kain-Fritsch积云对流参数化、WSM6 微物理方案、YSU 边界层方案、Dudhia 短波辐射、RRTM 长波辐射方案。背景场与边界条件取自NCEP 0.5°×0.5°分辨率的GFS全球预报场，同化了常规及加密地面观测资料，常规探空资料，船舶、浮标、测风、飞机、辽宁省自动气象站观测资料。此外，还实现了中国北方地区60 部雷达的径向风和反射率因子、葵花卫星云导风资料的同化。系统每日08、20 时（北京时，下同）启动两次，预报时效为84 h。

3 GNSS-ZTD(ZTD)资料3DVar 方案
3.1 ZTD 资料3DVar 方案简介
3DVar 同化目标函数定义为（Parrish，et al，1992）

pagenumber_ebook=111,pagenumber_book=828
式中，x0为大气状态变量，xb为大气状态的先验信息（称作背景场），yobs是已知观测向量，H为观测算子，B为 背景场（xb）的误差协方差矩阵（东北区域中尺度数值预报系统中由NMC 方法（Parrish，et al，1992）估计得到），O为 观测资料（yobs）的误差协方差矩阵，F 为观测算子（H）的误差协方差矩阵。目标函数式（1）的极小值是通过数学家提供的逐步迭代极小化方法得到的，其解即为分析场。

式（1）也可简化为

pagenumber_ebook=111,pagenumber_book=828
式中，x 为模式状态向量，Jb 、Jo分别为背景场和观测目标函数。

同化ZTD 观测需要在式（2 ） Jo项中加入ZTD 的贡献量

pagenumber_ebook=111,pagenumber_book=828
式中，HZTD 为由模式变量计算ZTD 的观测算子（详见3.2 节），R为ZTD 观测误差协方差矩阵（详见3.3 节），其他变量同式（1）。

3.2 ZTD 观测算子
ZTD 观测算子用于计算给定站点位置（包括高度）和时间条件下的模式ZTD 对应量，ZTD 计算由干延迟（ZHD）和湿延迟（ZWD）计算两部分组成，WRFDA 先在模式网格计算ZTD

pagenumber_ebook=111,pagenumber_book=828
式中，pa 为地面气压，θ 为纬度，ha 为地形高度（Saastamoinen，1972）。

pagenumber_ebook=111,pagenumber_book=828
式中，Q 为比湿，p 为气压，T 为温度，H 为高度，N 为模式面垂直层数，k′=2.21×10-7 K/Pa，k3=3.7×10-3 K2/Pa。

利用式（4）和（5）计算的模式格点ZTD（ZHD与ZWD 之和），先线性插值得到测站经纬度位置、模式地形高度处的ZTD 模式对应量，再由模式地形高度与测站高度的差异计算地形差补偿值对插值结果进行订正（仲跻芹等，2017），进而得到测站位置（包括高度）的模式ZTD 对应量。WRFDA 未考虑模式顶以上的大气延迟及大气折射率随高度减小的情况，这些可通过质量控制流程中偏差订正部分进行处理（Poli，et al，2007；Bennit，et al，2012）。

3.3 观测误差确定方法
3.3.1 传统方法

传统方法是根据华北区域ZTD 同化时采用的观测误差（15 mm）和2016 年6—8 月华北—辽宁区域经质量控制后所有站点的单站观测背景差月序列（以下称Somb）标准差最大值（26.5 mm），分别取15 和26.5 mm 观测误差进行对比试验。

3.3.2 Des 方法

Desroziers 等（2005）推导得出，在满足观测、背景、Omb、Oma 无偏、观测空间由估计的背景场误差协方差和观测误差协方差计算的增益矩阵与由真正的误差协方差得到的增益矩阵一致等条件下：

pagenumber_ebook=112,pagenumber_book=829
式中，E 为统计期望线性算子，pagenumber_ebook=112,pagenumber_book=829为新息增量矢量，pagenumber_ebook=112,pagenumber_book=829为观测分析差矢量，R为观测误差协方差。

采用式（7）（Desroziers，et al，2005）一次诊断观测误差，原因如前文所述。此外，考虑到文中使用的ZTD 观测覆盖范围广、地理环境差异较大、布设部门和仪器不统一、资料解算部门也不统一，导致站点间观测误差差异较为明显（详见图1），因此，参考Yan 等（2009a）的方法逐站进行观测误差统计。

pagenumber_ebook=112,pagenumber_book=829
式中，i 为观测站点序号，pagenumber_ebook=112,pagenumber_book=829为第i 个观测站点的观测误差诊断值，Si为i 站的样本数（本试验不小于100），pagenumber_ebook=112,pagenumber_book=829分别表示样本j 的观测值、与之对应的分析值和背景值，以下称之为Des 方法。

采用的Des 诊断流程：首先，将6—8 月Somb 标准差最大值26.5 mm 设为ZTD 观测误差假定值；其次，运行2016 年6—8 月6 h 间隔的WRF 3DVar 得到批量Omb 和Oma；再次，将批量Omb 和Oma 代入式（7）逐站诊断观测误差；最后，将诊断值写入对应的ZTD 观测文件，即可进行ZTD 同化与预报。

从华北—辽宁区域质量控制后保留的GNSS站点观测误差诊断值空间分布（图1）可知：观测误差诊断值明显小于假定值，假定值26.5 mm 时，诊断平均值为8.7 mm，诊断值与相同季节Bennitt 等（2012）、Poli 等（2007）确定的观测误差值（8 和10 mm）更为接近；观测站点间的诊断值差别较大，说明逐站进行观测误差诊断的必要性。此外，辽宁地区的观测误差诊断值总体小于华北地区，如4.1 节所述，这两部分资料来源不同，因此，该现象应源于资料解算与处理过程中的某些差异。本研究在采用较优的WRFDA 同化分析系统、基于NMC 方法（Parrish，et al，1992）估计背景场误差、由Omb 的统计确定观测误差假定值等条件下，理论上能够得到观测误差的合理估计（Ménard，2016），但由于无法获得观测真值，因此，无法确定诊断结果接近真值的程度。

pagenumber_ebook=113,pagenumber_book=830
图1 各GNSS 站点基于Des 方法的观测误差诊断值的空间分布
Fig.1 Spatial distribution of Des diagnosed observation errors for individual stations

4 GNSS-ZTD 资料及其在东北区域中尺度数值预报系统中的同化试验
4.1 ZTD 资料简介
使用的ZTD 资料来自京津冀晋鲁辽地区共279 个观测站，时间分辨率为1 h，华北地区ZTD 资料来源于国家气象信息中心采用GAMIT（King，et al，2000）软件解算及下发，辽宁省ZTD 资料是使用GAMIT 软件自行解算，站点高程（WGS84 高程）与东北区域中尺度数值预报系统地形高程一致。

根据中国辽宁和华北地区ZTD 资料实际情况，并结合以往ZTD 质量控制经验（Poli，et al，2007；Yan，et al，2009b；Bennitt，et al，2012；仲跻芹等，2017），在华北地区ZTD 质量控制流程（仲跻芹等，2017）基础上增加了高度差检查（如果模式面与GNSS 接收机高度差大于300 m，则剔除该站点）和空间稀疏化（阈值取25 km）步骤。

4.2 试验设计
已有研究（Macpherson，et al，2008；Vedel，et al，2004）表明，GNSS-ZTD 资料同化仅对可降水量大值区或降水过程影响明显，因此，选取2016 年汛期（6—8 月）华北—辽宁区域内13 d 的强降水个例，在东北区域中尺度数值预报系统平台上进行两种观测误差确定方法的ZTD 资料同化预报对比试验，本研究与实际业务运行不同的是，背景场与侧边界条件采用1°×1°的美国环境预报中心（NCEP）再分析资料，除全球卫星导航系统ZTD 外仅同化了常规地面观测资料，常规探空资料，船舶、浮标、测风、飞机探测资料，起报时间为08 时，积分24 h。同化对比试验除不同化ZTD的对照试验（Ctl）外，共设计3 个ZTD 观测误差确定方案进行对比试验：第一、二组ZTD 观测误差分别取15 mm （15 mm 方案）和26.5 mm （26.5 mm 方案），第三组试验是基于Des 方法诊断观测误差的Des 方案。

4.3 试验结果
4.3.1 同化分析效果检验

从质量控制流程中空间稀疏化步骤剔除的站点中选取了辽宁省12 个站点进行同化分析效果的检验，从对照试验及3 个同化方案的分析场模拟ZTD 相对于观测的平均绝对误差（图2）可知，13 个个例均为同化后误差小于未同化，除图中黑色虚框标识的4 个个例外，其他9 个个例（69%）均为Des方案误差最小。此外，平均误差分布（图略）也均为同化后平均误差更接近于0，Des 方案平均误差绝对值小于其他方案情况居多。表明ZTD 同化优化了分析场，Des 方法诊断的ZTD 观测误差更为合理、提高了分析场的精度。

4.3.2 降水预报效果评估

本研究基于模式二层嵌套范围内共计709 个常规降水观测，进行了08 时起报的24 h 降水预报ETS（Equitable Threat Score） 评分（Mesinger，2008）和预报偏差（预报站点数与实况站点数的比值）统计。雨量等级划分：小雨（及以上）量级：≥0.1 mm；中雨（及以上）量级：≥10 mm；大雨（及以上）量级：≥25 mm；暴雨（及以上）量级：≥50 mm；大暴雨（及以上）量级:≥100 mm。

pagenumber_ebook=114,pagenumber_book=831
图2 各试验方案ZTD 分析误差平均绝对值对比 （虚框标识Des 方案未体现优势）
Fig.2 Mean absolute errors of ZTD analysis（The bars surrounded by dashed lines indicate the situation that Des scheme didn't show advantage）

从ETS 评分和预报偏差检验结果（图3）可知，同化ZTD 资料后各降水量级（暴雨除外）的ETS 略有增大（图3a），中雨、大雨、大暴雨（以上）量级增加明显，对照试验小雨、中雨、大雨、暴雨、大暴雨（以上）量级的ETS 分别为0.31、0.29、0.23、0.18、0.10，Des 方案ETS 评分大多数高于传统方案，Des 方案与对照试验在各降水量级的ETS 评分差异分别为0.0037、0.0082、0.018、-0.005、0.018。从预报偏差（图3b）看，同化ZTD 资料后大暴雨（以上）量级降水预报强度有效增强，对应同化ZTD 后更高的ETS 评分，显示华北—辽宁地区ZTD 资料的同化可以提高模式对该区域降水尤其是强降水的强度、落区预报性能，Des 方案表现较优。

4.3.3 高空气象要素预报效果评估

本研究取模式二层嵌套范围内北京（54511）、张家口（54401）、乐亭（54539）、沈阳（54342）、大连（54662）5 个常规探空站500、700、850 hPa 相对湿度、温度、风速预报相对实况的13 d 平均绝对误差、平均相对误差进行不同方案的预报效果对比。

在相对湿度检验中（图4a、b），对照试验和各试验方案的平均相对误差基本均为负值，即对湿度的预报均偏低，同化ZTD 后，负偏差绝对值小于对照试验，且各预报时次的平均绝对误差均小于对照试验，其中Des 方案误差最小，表明同化ZTD 资料后，相对湿度预报与实况更接近，可信度更高；在温度的检验中（图4c、d），同化ZTD 资料后，略微改善了对照试验温度预报偏高或偏低的情况，且平均绝对误差基本小于对照试验，Des 方案误差也最小；在风速的检验中（图4e、f），对照试验500 hPa 的风速预报略偏小，而700、850 hPa 风速预报偏大，同化ZTD 资料后，正、负偏差均有所减小，平均绝对误差也小于对照试验，Des 方案误差也小于传统方案；相对湿度、温度、风速的平均绝对误差和平均相对误差的95%的置信区间（图4 中误差棒）均随预报时效的延长而略有增大，各要素同化ZTD 后95%置信区间较对照试验均略有减小。综上所述，同化GNSS-ZTD 资料后，湿度预报偏低等情况得到了一定的改善，提高了模式对高空常规要素的预报性能，Des 方案预报效果相对较好。

5 2016 年7 月25 日的个例分析
在中心位于贝加尔湖东侧的切断低压东移南下与副热带高压（副高）西进北上的共同作用下（图5a），2016 年7 月25 日凌晨开始东北地区中南部发生了一次大范围强降水过程，降水主要集中在25 日08 时—26 日08 时，100 mm 以上强降水区位于辽宁省东北部（图5b），降水中心抚顺市经济开发区站（下称抚顺开发区站）24 h 累计雨量为195.5 mm，副高的水汽输送向来是辽宁乃至东北地区区域性强降水过程发生的必要条件，26 日08 时后随着副高的南退及冷涡的东移强降水过程结束。

5.1 ZTD 3DVar 对初始场的影响
由25 日08 时—26 日08 时累计降水实况与各方案预报对比可知，Des 方案的预报与实况最为接近，能够修正对照试验在辽宁东部强雨区预报偏弱的情况，因此，以Des 方案为例，考察ZTD 资料同化对模式分析与预报的影响。

pagenumber_ebook=114,pagenumber_book=831
图3 13 d 平均的24 h 降水预报ETS 评分 （a） 与预报偏差 （b）
Fig.3 Averaged ETS （a） and bias （b） of 24 h precipitation forecasts for all the thirteen days

pagenumber_ebook=115,pagenumber_book=832
图4 3 种试验方案与对照试验的12 h （a、c、e） 和24 h （b、d、f） 相对湿度 （a、b）、温度 （c、d）、风速 （e、f） 预报平均绝对误差 （abias）、平均相对误差 （bias）（c—f 左侧坐标为平均绝对误差，右侧坐标为平均相对误差，误差棒为平均绝对误差和平均相对误差的95%置信区间）
Fig.4 12 h （a，c，e） and 24 h （b，d，f） forecast mean absolute biases and mean biases of relative humidity （a，b），temperature （c，d），wind speed （e，f） from three test schemes and the Ctl experiment verified against observations（the left and right axes in Figure c—f are for mean absolute bias and bias respectively；the error bars show the 95%confidence intervals of average absolute bias and average bias）

首先考察初始时刻Des 方案与对照试验 700 hPa的比湿差值分布（图6a），可见，ZTD 观测区域（图1）内除辽宁中部出现狭长弱负值区外，其余基本均为较大的正差值区，最大增量达0.0015 g/g 以上，位于辽宁南部和河北西北部（图6a 中三角形与圆圈标记处），其他对流层中低层比湿场也有类似改变。对比ERA5 再分析与对照试验的700 hPa 比湿差值（图6b）可知，两比湿差值场数值分布相似，均在辽宁南部、河北北部出现正差值区，正差值中心（图6b中三角形与圆圈标记处）量级相当，在河北南部、辽宁西南部均出现负差值区，但范围、强度有差别，其他对流层中低层的两比湿差值场也对应较好。说明Des 方案对对照试验进行了有效修正，同化ZTD资料后雨区及其上游分析场获得了有效增湿。

pagenumber_ebook=116,pagenumber_book=833
图5 2016 年7 月25 日 （a） 14 时500 hPa 高度 （等值线，单位：dagpm）、700 hPa 风和850 hPa 相对湿度 （阴影，单位：%） 场，（b） 25 日08 时—26 日08 时累计降水量 （单位：mm，星形标记降水中心抚顺开发区站）
Fig.5 （a） 500 hPa height （contours，unit:dagpm），700 hPa wind and 850 hPa relative humidity （shadings，unit:%）at 14:00 BT 25 July 2016；（b） observational precipitation from 08:00 BT 25 July to 08:00 BT 26 July 2016（unit:mm，star marks precipitation center in Fushun development zone station）

pagenumber_ebook=116,pagenumber_book=833
图6 2016 年7 月25 日08 时700 hPa 比湿分析场差值 （单位：g/g） 分布（a.Des 方案与对照试验差值（黑色虚线为后续分析剖面），b.ERA5 再分析与对照试验差值；圆圈、三角形标记两个正差值中心）
Fig.6 Differences of specific humidity （unit:g/g） analysis in 700 hPa at 08:00 BT 25 July 2016 （a.Des-Ctl （the dashed black line is where the cross-section is made），b.ERA5-Ctl；the circle and triangle represent the centers of positive difference）

pagenumber_ebook=116,pagenumber_book=833
图7 2016 年7 月25 日08 时比湿 （单位：g/g） 垂直分布 （a.模式格点（40°N，123°E）比湿廓线分布（黑色实线：ERA5 再分析，红色虚线：Des 方案，蓝色虚线：对照试验），b.沿图6a 中黑色虚线的比湿差值垂直剖面（阴影为ERA5 再分析与对照试验的差值为正的部分，等值线为Des 方案与对照试验的差值））
Fig.7 Vertical distribution of specific humidity （unit:g/g） at 08:00 BT 25 July 2016 （a.specific humidity profile at （40°N，123°E），black solid line:ERA5 reanalysis，red dotted line:Des scheme，blue dotted line:Ctl experiment；b.cross-section of specific humidity difference along the black dashed line shown in Fig.6a，shadings show above-zero differences between ERA5 reanalysis and Ctl experiment，and contours indicate the differences between Des results of scheme and Ctl experiment）

其次，考察初始时刻700 hPa 比湿增值中心处（40°N，123°E，即图6a 中三角形标记处）Des 方案与对照试验比湿垂直廓线分布对比，结果表明，在该处400 hPa 以下Des 方案比湿分析廓线较对照试验明显偏湿（图7a）、更接近ERA5 再分析廓线。此外，从经过该处沿图6a 中黑色虚线的垂直剖面上两方案的比湿分布对比（图7b）可知，该剖面上ERA5再分析减对照试验比湿主要为正差值，大差值区位于122.6°—124°E 的900—650 hPa 高度，而该大差值区内，Des 方案减对照试验也为正值，且两正差值中心基本重合、数值相近。总之，同化ZTD 资料后，不仅700 hPa 比湿增值中心处中低层比湿廓线有效增湿，经过该处中低层水汽输送通道沿线湿度也获得了有效提高。

温、压、风分析场在同化前后未发生明显改变，ZTD 资料同化主要影响的是湿度场，东北区域中尺度数值预报系统3DVar 的控制变量选项CV5中湿度与其他控制变量均不相关，这应该是造成同化后温、压、风分析场看不出明显变化的主要原因。

5.2 ZTD 3DVar 对积分初期水成物分布及降水模拟效果的影响
以第2 小时（10 时）云水、云冰、雨水、雪、霰含量沿图6a 黑色虚线的垂直分布作为代表进行分析，图8 给出对照试验和Des 方案对该时次各类水成物混合比的模拟（图8b、d、f—h），同时给出ERA5的云水、云冰、雨水含量（可近似表示混合比）再分析（图8a、c、e）作为参考。由图8a 可知，ERA5 云水含量主要分布于123.3°—124.8°E 的 800—1000 hPa、中心极值约0.0004 g/g，Des 方案主要分布位置、中心极值均与再分析接近（图8b 阴影），而对照试验比再分析偏西南、分布范围小（图8b 等值线）；ERA5 云冰含量主要分布于123.5°—124.8°E 的250—500 hPa、中心极值约0.00012 g/g（图8c），Des 方案的位置、范围及强度均较对照试验更为接近再分析（图8d）；ERA5 雨水含量较低、中心极值约0.00003 g/g（图8e），从分布位置看，Des 方案更接近再分析（图8f）；对照试验云水、雨水分布位置上空无雪、霰分布，而Des 方案雪、霰均有较大范围的分布，中心极值分别约为0.00045 和0.0035 g/g（图8g、h），无ERA5 再分析资料进行对比。积分初期其他时次该剖面及附近水成物分布也均为Des 方案较对照试验含量、范围偏大，位置偏东北、更接近再分析。

综上所述，同化ZTD 资料后，副热带高压水汽输送通道沿线降水有关水成物的含量有效提高、分布位置向东北方向偏移约1 个纬度，这势必带来降雨量和雨区位置预报的调整。由2016 年7 月25 日08—12 时地面降水分布实况（图9a1—a4）可见，降水主要分布在中国蒙冀辽交界地带及辽宁东部地区。对比两方案预报可知，对照试验和Des 方案降水预报差别主要发生在辽宁东部地区：在第1 小时，对照试验未预报出辽宁东部降水（图9b1），而Des 方案报出零散的降水（图9c1），在随后的模拟中，对照试验仅在第3—4 小时预报辽宁东部有小范围弱降水（图9b3—b4），而Des 方案在第2 小时即开始预报辽宁东部明显的强降水（图9c2），预报范围和强度在第3—4 小时明显增大、增强，至11—12 时已接近实况（图9c4）。

5.3 ZTD 3DVar 对积分初期温度模拟效果的影响
降水过程通过凝结潜热释放加热大气进而影响气温的分布（丁一汇，2005），因此，对照试验和Des 方案降水模拟的差异会带来温度预报效果的不同。从第2 个积分时次对照试验（等值线）与Des 方案（阴影）沿图6a 黑色虚线的凝结潜热释放垂直剖面（图10a）可知，潜热释放位置与降水位置对应，两方案大值中心分别位于122.3°E 和123.7°E 附近，Des 方案潜热释放强度和范围均大于对照试验。两方案均在750 hPa 高度（图10a 中虚线）附近潜热释放强度最大，从该高度处的温度预报与再分析对比（图10b—d）可知，Des 方案有效地向东北方向推进了286—288 K 等温线，在辽宁东南部123°E 附近，对照试验287 K 等温线紧贴海岸线（图10c），而Des 方案（图10d）该等温线已深入内陆、更接近再分析。类似的温度模拟效果改进也发生在850—600 hPa 的其他高度和积分初期的其他时次。

5.4 ZTD 3DVar 对积分初期风场模拟效果的影响
降水过程也通过潜热释放激发地转偏差（实际风与地转风的矢量差）影响降水区附近的流场（李毓芳，1982；倪允琪等，2006），因此，对照试验和Des 方案降水模拟的差异也会带来风场预报效果的不同。首先，由图11a、b 可知，Des 方案与对照试验地转偏差流型相似，地转偏差风速（图11a 阴影）大于5 m/s 的区域与降水落区基本对应，说明该地转偏差主要由降水过程中的潜热释放所激发；其次，Des 方案地转偏差风速在河北、辽宁的大部分地区均大于对照试验（图11b 阴影），导致两方案在该处的预报风速差大值区（图11c）；最后，由两方案预报风场与ERA5 再分析的对比（图11d—f）可知，Des 方案与对照试验17 m/s 以上风速大值区均较ERA5 明显偏西北、偏强（图11d—f），但相比之下，Des 方案该风速大值区（在全球卫星导航系统ZTD 同化区域内的）东、南侧边界较对照试验有明显的向东、向南扩大、向再分析场靠近的趋势，Des 方案风场模拟较对照试验向再分析场靠近的趋势也发生在对流层中低层的其他高度和积分初期的其他时次。

pagenumber_ebook=119,pagenumber_book=836
图8 2016 年7 月25 日10 时沿图6a 黑色虚线的水成物 （单位：g/g） 分布垂直剖面（a、c、e.ERA5 再分析云水、云冰、雨水，b、d、f、g、h.对照试验 （等值线） 与Des 方案 （色阶） 模拟的云水、云冰、雨水、雪、霰）
Fig.8 Cross-sections of hydrometeors （unit:g/g） distribution along the dashed line shown in Fig.6a at 10:00 BT 25 July 2016 （a，c，e.the ERA5 reanalysis of cloud water，cloud ice，rain ；b，d，f，g，h.cloud water，cloud ice，rain，snow，graupel simulated by Ctl experiment （contours） and Des scheme （shadings））

pagenumber_ebook=119,pagenumber_book=836
图9 实况 （a1—a4）、对照试验 （b1—b4） 和Des 方案 （c1—c4） 积分前4 h （2016 年7 月25 日08—12 时） 逐时累计降水量（a1—c1.08：00—09：00 BT，a2—c2.08：00—10：00 BT，a3—c3.08：00—11：00 BT，a4—c4.08：00—12：00 BT；单位：mm）
Fig.9 Accumulated rainfall （unit:mm） of observation and simulation of Ctl experiment and Des scheme at every hour over the period from 08:00 to 12:00 BT 25 July 2016 （a1—c1.08：00—09：00 BT，a2—c2.08：00—10：00 BT，a3—c3.08：00—11：00 BT，a4—c4.08：00—12：00 BT）

6 结论与讨论
基于2016 年6—8 月华北—东北地区的地基GNSS-ZTD 观测资料、东北区域中尺度数值预报模式及其三维变分同化系统，以2016 年6—8 月13 d的强降水为例，开展了基于Des 方法和传统方法确定观测误差的ZTD 资料同化预报批量对比试验，并对分析、预报结果进行了统计检验和典型个例分析，结果表明：

（1）基于本研究确定的Des 方案诊断流程得到的ZTD 观测误差诊断值较为合理，诊断值站点间差异较大，说明逐站诊断观测误差的必要性。

（2）批量分析及预报效果检验表明，华北—辽宁地区ZTD 资料的同化可以使模式对该区域降水尤其是中雨以上量级降水的强度、落区预报性能获得提高；使温、湿、风等要素的预报向观测靠近，Des 方案分析、预报效果优于传统方案。

（3）2016 年7 月25 日ZTD 同化预报效果分析表明，ZTD 资料同化有效增强了初始湿度场，对温、压、风初值无明显的调整，改进了积分初期水成物含量与分布位置进而改善了降水预报效果，修正了对照试验对辽宁东部地区强降水的明显漏报，且通过降水的反馈作用改进了温度与风场预报。

pagenumber_ebook=121,pagenumber_book=838
图10 积分第2 小时 （2016 年7 月25 日10 时）（a）对照试验 （等值线） 与Des 方案 （色阶） 沿图6a 虚线的凝结潜热 （单位：J/（g·h），虚线标识750 hPa 高度） 垂直剖面和 （b） ERA5 再分析、（c） 对照试验与（d） Des 方案750 hPa 温度 （单位：K） 分布 （虚线标识286、287、288 K 等温线）
Fig.10 Cross-sections of condensational latent heat at the second hour of integration （10:00 BT 25 July 2016） from （a） Ctl（contours） and Des （shadings） along the dashed line shown in Fig.6a （unit:J/（g·h）；dashed line marks the height of 750 hPa） and temperature at 750 hPa from （b） ERA5，（c） Ctl，and （d） Des （unit:K，dashed lines mark the isotherms of 286，287，288 K）

pagenumber_ebook=121,pagenumber_book=838
图11 积分第3 小时 （2016 年7 月25 日11 时） 700 hPa 风场 （色阶，单位：m/s；a.Des 方案地转偏差流场与地转偏差风速，b.对照试验地转偏差流场与Des-Ctl 地转偏差风速，c.Des-Ctl 风速，d.ERA5 再分析风场，e.对照试验风场，f.Des 方案风场）
Fig.11 700 hPa wind field （shadings，unit:m/s） at the third hour of integration at 11:00 BT 25 July 2016（a.streamline and speed of geostrophic wind deviations of Des，b.streamline of geostrophic wind deviation of Ctl and speed difference of geostrophic wind deviation between Des and Ctl，c.wind speed difference between Des and Ctl，d.ERA5 wind field，e.Ctl wind field，f.Des wind field）

观测误差的确定是一项较为复杂的工作，已有文献（Daescu，et al，2010，2013）表明，模式预报对观测误差协方差存在敏感性，提出将诊断与（根据敏感性）调整相结合是更为理想的观测误差确定方式，这也是值得继续深入研究的问题。

参考文献

陈敏,范水勇,仲跻芹等.2010.全球定位系统的可降水量资料在北京地区快速更新循环系统中的同化试验.气象学报，68（4）：450-463.Chen M,Fan S Y,Zhong J Q,et al.2010.An experimental study of assimilating the Global Position System-precipitable water vapor observations into the rapid updated cycle system for the Beijing area.Acta Meteor Sinica，68（4）：450-463 （in Chinese）

丁金才.2009.GPS 气象学及其应用.北京:气象出版社,265pp.Ding J C.2009.GPS Meteorology and its Application.Beijing:China Meteorological Press,265pp (in Chinese)

丁一汇.2005.高等天气学.第2 版.北京:气象出版社,585pp.Ding Y H.2005.Advanced Synoptic Meteorology.2nd ed.Beijing: China Meteorological Press,585pp (in Chinese)

李毓芳.1982.暴雨过程中的三支气流:一次暴雨过程的初步分析.杭州大学学报，9（2）：207-218.Li Y F.1982.The three air streams:The SESAME case Ⅲ,1979.J Hangzhou Univ，9（2）：207-218

倪允琪,周秀骥.2006.“我国重大天气灾害形成机理与预测理论研究”取得的主要研究成果.地球科学进展，21（9）：881-894.Ni Y Q,Zhou X J.2006.Main scientific issues and achievements of state 973 project on study for formation mechanism and prediction theories of severe weather disasters in China.Adv Earth Sci，21（9）：881-894 （in Chinese）

王金成,龚建东,赵滨.2015.一种新的COSMIC 大气折射率资料观测误差估计方法及在GRAPES 全球三维变分同化中的应用.气象学报，73（1）：142-158.Wang J C,Gong J D,Zhao B.2015.A new method for estimating observation error of the COSMIC refractivity data and its impacts on GRAPES-GFS model weather forecasts.Acta Meteor Sinica，73（1）：142-158 （in Chinese）

仲跻芹,GUO Y R,张京江.2017.华北地区地基GPS 天顶总延迟观测的质量控制和同化应用研究.气象学报，75（1）：147-164.Zhong J Q,Guo Y R,Zhang J J.2017.A study of quality control and assimilation of groundbased GPS ZTD in North China.Acta Meteor Sinica，75（1）：147-164

Arriola J S,Lindskog M,Thorsteinsson S,et al.2016.Variational bias correction of GNSS ZTD in the HARMONIE modeling system.J Appl Meteor Climatol，55（5）：1259-1276

Bennitt G V,Jupp A.2012.Operational assimilation of GPS zenith total delay observations into the Met Office numerical weather prediction models.Mon Wea Rev，140（8）：2706-2719

Boniface K,Ducrocq V,Jaubert G,et al.2009.Impact of high-resolution data assimilation of GPS zenith delay on Mediterranean heavy rainfall forecasting.Ann Geophys，27（7）：2739-2753

Chapnik B,Desroziers G,Rabier F,et al.2004.Properties and first application of an error-statistics tuning method in variational assimilation.Quart J Roy Meteor Soc，130（601）：2253-2275

Chapnik B,Desroziers G,Rabier F,et al.2006.Diagnosis and tuning of observational error in a quasi-operational data assimilation setting.Quart J Roy Meteor Soc，132（615）：543-565

Daescu D N,Todling R.2010.Adjoint sensitivity of the model forecast to data assimilation system error covariance parameters.Quart J Roy Meteor Soc，136（653）：2000-2012

Daescu D N,Langland R H.2013.Error covariance sensitivity and impact estimation with adjoint 4D-Var:Theoretical aspects and first applications to NAVDAS-AR.Quart J Roy Meteor Soc，139（670）：226-241

Dee D P.1995.On-line estimation of error covariance parameters for atmospheric data assimilation.Mon Wea Rev，123（4）：1128-1145

Dee D P,Da Silva A M.1999.Maximum-likelihood estimation of forecast and observation error covariance parameters.Part I:Methodology.Mon Wea Rev，127（8）：1822-1834

Desroziers G,Ivanov S.2001.Diagnosis and adaptive tuning of observationerror parameters in a variational assimilation.Quart J Roy Meteor Soc，127（574）：1433-1452

Desroziers G,Berre L,Chapnik B,et al.2005.Diagnosis of observation,background and analysis-error statistics in observation space.Quart J Roy Meteor Soc，131（613）：3385-3396

Giannaros C,Kotroni V,Lagouvardos K,et al.2020.Assessing the impact of GNSS ZTD data assimilation into the WRF modeling system during high-impact rainfall events over GREECE.Remote Sens，12（3）：383

Hollingsworth A,Lönnberg P.1986.The statistical structure of short-range forecast errors as determined from radiosonde data.Part I:The wind field.Tellus A，38（2）：111-136

King B R W,Bock Y.2000.Documentation for the GAMIT GPS Analysis Software.Cambridge:Massachusetts Institute of Technology,206pp

Kuo Y H,Wee T K,Sokolovskiy S,et al.2004.Inversion and error estimation of GPS radio occultation data.J Meteor Soc Jpn，82（1B）：507-531

Lindskog M,Ridal M,Thorsteinsson S,et al.2017.Data assimilation of GNSS zenith total delays from a Nordic processing centre.Atmos Chem Phys，17（22）：13983-13998

Macpherson S R,Deblonde G,Aparicio J M,et al.2008.Impact of NOAA ground-based GPS observations on the Canadian regional analysis and forecast system.Mon Wea Rev，136（7）：2727-2746

Mascitelli A,Federico S,Fortunato M,et al.2019.Data assimilation of GPSZTD into the RAMS model through 3D-Var:Preliminary results at the regional scale.Meas Sci Technol,doi:10.1088/1361-6501/ab0b87

Ménard R.2016.Error covariance estimation methods based on analysis residuals:Theoretical foundation and convergence properties derived from simplified observation networks.Quart J Roy Meteor Soc，142（694）：257-273

Mesinger F.2008.Bias adjusted precipitation threat scores.Adv Geosci，16：137-142

Parrish D F,Derber J C.1992.The National Meteorological Center's spectral statistical-interpolation analysis system.Mon Wea Rev，120（8）：1747-1763

Poli P,Moll P,Rabier F,et al.2007.Forecast impact studies of zenith total delay data from European near real-time GPS stations in Météo France 4DVAR.J Geophys Res，112（D6）：D06114

Rohm W,Guzikowski J,Wilgan K,et al.2019.4DVAR assimilation of GNSS zenith path delays and precipitable water into a numerical weather prediction model WRF.Atmos Meas Tech，12（1）：345-361

Saastamoinen J.1972.Atmospheric correction for the troposphere and stratosphere in radio ranging satellites∥Henriksen S W,Mancini A,Chovitz B H.The Use of Artificial Satellites for Geodesy.Washington:American Geophysical Union,247-251

Todling R.2015.A complementary note to "a lag-1 smoother approach to system-error estimation":The intrinsic limitations of residual diagnostics.Quart J Roy Meteor Soc，141（692）：2917-2922

Vedel H,Mogensen K S,Huang X Y.2001.Calculation of zenith delays from meteorological data comparison of NWP model,radiosonde and GPS delays.Phys Chem Earth A Solid Earth Geodesy，26（6-8）：497-502

Vedel H,Huang X Y.2004.Impact of ground based GPS data on numerical weather prediction.J Meteor Soc Jpn，82（1B）：459-472

Waller J A,Dance S L,Nichols N K.2016.Theoretical insight into diagnosing observation error correlations using observation-minus-background and observation-minus-analysis statistics.Quart J Roy Meteor Soc，142（694）：418-431

Yan X,Ducrocq V,Jaubert G,et al.2009a.The benefit of GPS zenith delay assimilation to high-resolution quantitative precipitation forecasts:A case study from COPS IOP 9.Quart J Roy Meteor Soc，135（644）：1788-1800

Yan X,Ducrocq V,Poli P,et al.2009b.Impact of GPS zenith delay assimilation on convective-scale prediction of Mediterranean heavy rainfall.J Geophys Res，114（D3）：D03104

Yang S C,Huang Z M,Huang C Y,et al.2020.A case study on the impact of ensemble data assimilation with GNSS-zenith total delay and radar data on heavy rainfall prediction.Mon Wea Rev，148（3）：1075-1098

GONG Ying1 LI Deqin1 ZHONG Jiqin2 ZHOU Xiaoshan1 YANG Sen1 PAN Xiao1

1.Institute of Atmosphere Environment， CMA， Shenyang 110166， China
2.Institute of Urban Meteorology， CMA， Beijing 100089， China

Abstract The Regional Meso-scale Numerical Prediction System of Northeast China (RMNPSNC) is used to study the application of an observation error diagnostic method based on Desroziers et al (2005) (here in after Des method) in the GNSS-ZTD (Global Navigation Satellite System-Zenith Total Delay,ZTD) 3DVar assimilation.The Des method is compared with traditional observation error determination method based on ZTD assimilation and forecast experiments of thirteen rainfall cases during the period from June to August 2016.A comparative study is then conducted based on assimilation and forecast results with and without ZTD data to evaluate the effect of assimilating ZTD data into the RMNPSNC.The results are summarized as follows.(1) The errors of ZTD observations diagnosed by the Des method are relatively reasonable,and the differences of diagnostic values between stations are large,indicating the necessity of diagnosing observation errors station by station.(2) Assimilation of the ZTD data improves the forecast of the intensity and distribution of heavy rain and the forecasts of temperature,humidity and wind are closer to observations.The analysis and forecast effects of the Des scheme are better than that of the traditional scheme.(3) For the heavy rain process in Northeast China on 25 July 2016,the assimilation of ZTD effectively increases the initial humidity field,improves the content and spatial distribution of hydrometeors at the initial hours of integration,corrects the failed precipitation forecast in the east of Liaoning province,and improves temperature and wind forecasts due to realistic precipitation feedback.The diagnostic ZTD observation errors obtained by the Des method are more reasonable than that by traditional method.Therefore,the assimilation and forecast can be improved by using the Des method.Progress in the forecast of rain,temperature,humidity and wind can be made through assimilating the ZTD.

Key words Observation error determination，Ground-based GNSS-ZTD，Data assimilation，3DVar

2019-07-16 收稿，2020-05-05 改回.

公颖，李得勤，仲跻芹，周晓珊，杨森，潘晓.2020.不同观测误差确定方法对地基GNSS-ZTD 资料同化预报效果影响的对比分析.气象学报，78（5）：826-839

Gong Ying，Li Deqin，Zhong Jiqin，Zhou Xiaoshan，Yang Sen，Pan Xiao.2020.A comparative analysis of influences of two different observation error determination methods on the assimilation of ground-based GNSS-ZTD data.Acta Meteorologica Sinica，78（5）:826-839

*资助课题：国家重点研发计划项目（2018YFC1506803、2017YFC1501806）、辽宁省气象局科学技术研究项目（201703）、公益性行业（气象）科研专项（GYHY201406012）、国家自然科学基金项目（41705011）。

作者简介：公颖，主要从事资料同化与中尺度数值模拟研究。E-mail：gongying74@sohu.com

通信作者：周晓珊，主要从事天气分析、资料同化等方面的研究。E-mail：xiaoshan_zhou@163.com

中图法分类号 P407