气象学报  2020, Vol. 78 Issue (5): 826-839   PDF    
http://dx.doi.org/10.11676/qxxb2020.053
中国气象学会主办。
0

文章信息

公颖, 李得勤, 仲跻芹, 周晓珊, 杨森, 潘晓. 2020.
GONG Ying, LI Deqin, ZHONG Jiqin, ZHOU Xiaoshan, YANG Sen, PAN Xiao. 2020.
不同观测误差确定方法对地基GNSS-ZTD资料同化预报效果影响的对比分析
A comparative analysis of influences of two different observation error determination methods on the assimilation of ground-based GNSS-ZTD data
气象学报, 78(5): 826-839.
Acta Meteorologica Sinica, 78(5): 826-839.
http://dx.doi.org/10.11676/qxxb2020.053

文章历史

2019-07-16 收稿
2020-05-05 改回
不同观测误差确定方法对地基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)数据    资料同化    三维变分    
A comparative analysis of influences of two different observation error determination methods on the assimilation of ground-based GNSS-ZTD data
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    
1 引 言

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

观测误差的确定是目前天顶总延迟资料同化工作中的关键与难点之一(Poli,et al,2007Macpherson,et al,2008Boniface,et al,2009Bennitt,et al,2012仲跻芹等,2017)。变分同化过程中,观测误差用于确定分析中观测资料的权重(Macpherson,et al,2008Waller,et al,2016),是观测信息能否被充分有效利用以及分析场能否对大气状态进行最优估计的重要决定因素,由于观测真值无法获得,观测误差的确定成为当前资料同化中一个主要的困难与挑战(Chapnik,et al,2004)。到目前为止,绝大多数的天顶总延迟资料同化中,观测误差的确定方法均较为粗略(以下称传统方法):通常基于天顶总延迟观测背景差(以下称Omb)标准差确定观测误差(Boniface,et al,2009Yan,et al,2009aBennitt,et al,2012);也有根据天顶总延迟平均值与Omb标准差间的稳定拟合关系确定观测误差(Macpherson,et al,2008);更多的是基于经验直接给定观测误差(Vedel,et al,2004Poli,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,19951999),该方法基于新息增量,通过最大似然估计确定误差协方差模式参数,进而得到观测误差协方差。(3)Desroziers等(2001)基于观测目标函数最小值的期望值与实际值来计算观测误差调整参数;此方法仅适用于变分分析框架。(4)Desroziers等(2005)提出基于Omb与观测分析差(以下简称Oma)诊断观测误差协方差的方法,详见3.3.2节;其优点是适用于任何分析框架。方法(4)是目前国际上较为流行的观测误差诊断方法,该方法应用于变分同化系统时需采用一次诊断的方式,而非理论上应采用的不动点迭代方法(Desroziers,et al,2005),因为在相关的观测误差被假定为不相关的变分同化系统中迭代方法不可行(Waller,et al,2016Todling,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

$\begin{split} { J}\left({{{ x}_0}} \right) =& \frac{\,1\,}{\,2\,}{\left({{{ x}_0} - {{ x}_{\rm{b}}}} \right)^{\rm{T}}}{{{B}}^{ - 1}}\left({{{ x}_0} - {{ x}_{\rm{b}}}} \right) +\\& \frac{\,1\,}{\,2\,}{\left({H\left({{{ x}_0}} \right) - {{ y}^{{\rm{obs}}}}} \right)^{\rm{T}}}{({{O}} + {{F}})^{ - 1}}\left({H\left({{{x}_0}} \right) - {{y}^{{\rm{obs}}}}} \right) \end{split}$ (1)

式中, ${{x}}_{0}$ 为大气状态变量, ${{x}}_{\rm{b}}$ 为大气状态的先验信息(称作背景场), ${{y}}^{\rm{obs}}$ 是已知观测向量, $ H $ 为观测算子, $ {{B}} $ 为背景场( ${{x}}_{\rm{b}}$ )的误差协方差矩阵(东北区域中尺度数值预报系统中由NMC方法(Parrish,et al,1992)估计得到), $ {{O}} $ 为观测资料( ${{y}}^{\rm{obs}}$ )的误差协方差矩阵, $ {{F}} $ 为观测算子( $ H $ )的误差协方差矩阵。目标函数式(1)的极小值是通过数学家提供的逐步迭代极小化方法得到的,其解即为分析场。

式(1)也可简化为

$ { J}\left({ x} \right) = {{ J}_{\rm{b}}}\left({ x} \right) + {{{ {J}}}_{\rm{o}}}\left({ x} \right) $ (2)

式中,x模式状态向量, ${{ J}}_{\rm{b}}$ ${{ J}}_{\rm{o}}$ 分别为背景场和观测目标函数。

同化ZTD观测需要在式(2) ${{ J}}_{\rm{o}}$ 项中加入ZTD的贡献量

$ { J}_0^{{\rm{ZTD}}}\left({{{ x}_0}} \right) = {\left({{H^{{\rm{ZTD}}}}\left({{{ x}_0}} \right) - {{ y}^{{\rm{obs}}}}} \right)^{\rm{T}}}{{{R}}^{ - 1}}\left({{H^{{\rm{ZTD}}}}\left({{{ x}_0}} \right) - {{ y}^{{\rm{obs}}}}} \right) $ (3)

式中, $ {H}^{\rm{ZTD}} $ 为由模式变量计算ZTD的观测算子(详见3.2节), $ {{R}} $ 为ZTD观测误差协方差矩阵(详见3.3节),其他变量同式(1)。

3.2 ZTD观测算子

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

$ {\rm{ZHD}} = \frac{{2.2768 \times {{10}^{ - 5}} \times {p_{\rm a}}}}{{1 - 2.66 \times {{10}^{ - 3}} \times{\rm{cos}} (2\theta) - 2.8 \times {{10}^{ - 7}} \times {h_{\rm a}}}} $ (4)

式中,pa为地面气压,θ为纬度,ha为地形高度(Saastamoinen,1972)。

$ {\rm{ZWD}} = \frac{1}{{0.622}}\sum\limits_{{{i}} = 1}^N {\frac{{{Q_{{i}}}{p_{{i}}}}}{{{T_{{i}}}}}} \left({{{k'}} + \frac{{{{{k}}_3}}}{{{T_{{i}}}}}} \right)\left({{H_{{{i}} + 1}} - {H_{{i}}}} \right) $ (5)

式中,Q为比湿,p为气压,T为温度,H为高度,N为模式面垂直层数, ${{k^{'}}} = 2.21 \times {10^{ - 7}}\;{\rm{K}}/{\rm{Pa}}$ ${{{k}}_3} = 3.7 \times {10^{ - 3}}\;{{\rm{K}}^2}/{\rm{Pa}}$

利用式(4)和(5)计算的模式格点ZTD(ZHD与ZWD之和),先线性插值得到测站经纬度位置、模式地形高度处的ZTD模式对应量,再由模式地形高度与测站高度的差异计算地形差补偿值对插值结果进行订正(仲跻芹等,2017),进而得到测站位置(包括高度)的模式ZTD对应量。WRFDA未考虑模式顶以上的大气延迟及大气折射率随高度减小的情况,这些可通过质量控制流程中偏差订正部分进行处理(Poli,et al,2007Bennit,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)推导得出,在满足观测、背景、OmbOma无偏、观测空间由估计的背景场误差协方差和观测误差协方差计算的增益矩阵与由真正的误差协方差得到的增益矩阵一致等条件下:

$ E\left[ {{{d}}_{\rm{a}}^{\rm{o}}{{\left({{{d}}_{\rm{b}}^{\rm{o}}} \right)}^{\rm{T}}}} \right] = {{R}} $ (6)

式中,E为统计期望线性算子, ${{d}}_{\rm{b}}^{\rm{o}}$ 为新息增量矢量, ${{d}}_{\rm{a}}^{\rm{o}}$ 为观测分析差矢量, $ {{R}} $ 为观测误差协方差。

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

图 1  各GNSS站点基于Des方法的观测误差诊断值的空间分布 Fig. 1  Spatial distribution of Des diagnosed observation errors for individual stations
$ {\widetilde {\left({\sigma _{{i}}^{\rm{o}}} \right)}^2} = \left({{{d}}_{\rm{a}}^{\rm{o}}} \right)_{{i}}^{\rm{T}}{\left({{{d}}_{\rm{b}}^{\rm{o}}} \right)_{{i}}}/{p_{{i}}} = \sum\limits_{{ j} = 1}^{{S_{{i}}}} {\left({y_j^{\rm{o}} - y_j^{\rm{a}}} \right)} \left({y_j^{\rm{o}} - y_j^{\rm b}} \right)/{p_{{i}}} $ (7)

式中,i为观测站点序号, $\widetilde {\left({\sigma _{{i}}^{\rm{o}}} \right)}$ 为第i个观测站点的观测误差诊断值, ${S_{{i}}}$ i站的样本数(本试验不小于100), $ {y}_{{j}}^{\rm{o}} $ $ {y}_{{j}}^{\rm{a}} $ $ {y}_{{j}}^{\rm{b}} $ 分别表示样本j的观测值、与之对应的分析值和背景值,以下称之为Des方法。

采用的Des诊断流程:首先,将6—8月Somb标准差最大值26.5 mm设为ZTD观测误差假定值;其次,运行2016年6—8月6 h间隔的WRF 3DVar得到批量OmbOma;再次,将批量OmbOma代入式(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),但由于无法获得观测真值,因此,无法确定诊断结果接近真值的程度。

4 GNSS-ZTD资料及其在东北区域中尺度数值预报系统中的同化试验 4.1 ZTD资料简介

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

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

4.2 试验设计

已有研究(Macpherson,et al,2008Vedel,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观测误差更为合理、提高了分析场的精度。

图 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)
4.3.2 降水预报效果评估

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

从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方案表现较优。

图 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
4.3.3 高空气象要素预报效果评估

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

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

图 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)
5 2016年7月25日的个例分析

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

图 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)
5.1 ZTD 3DVar对初始场的影响

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

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

图 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)

其次,考察初始时刻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比湿增值中心处中低层比湿廓线有效增湿,经过该处中低层水汽输送通道沿线湿度也获得了有效提高。

图 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)

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

5.2 ZTD 3DVar对积分初期水成物分布及降水模拟效果的影响

以第2小时(10时)云水、云冰、雨水、雪、霰含量沿图6a黑色虚线的垂直分布作为代表进行分析,图8给出对照试验和Des方案对该时次各类水成物混合比的模拟(图8bdfh),同时给出ERA5的云水、云冰、雨水含量(可近似表示混合比)再分析(图8ace)作为参考。由图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(图8gh),无ERA5再分析资料进行对比。积分初期其他时次该剖面及附近水成物分布也均为Des方案较对照试验含量、范围偏大,位置偏东北、更接近再分析。

图 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))

综上所述,同化ZTD资料后,副热带高压水汽输送通道沿线降水有关水成物的含量有效提高、分布位置向东北方向偏移约1个纬度,这势必带来降雨量和雨区位置预报的调整。

由2016年7月25日08—12时地面降水分布实况(图9a1a4)可见,降水主要分布在中国蒙冀辽交界地带及辽宁东部地区。对比两方案预报可知,对照试验和Des方案降水预报差别主要发生在辽宁东部地区:在第1小时,对照试验未预报出辽宁东部降水(图9b1),而Des方案报出零散的降水(图9c1),在随后的模拟中,对照试验仅在第3—4小时预报辽宁东部有小范围弱降水(图9b3b4),而Des方案在第2小时即开始预报辽宁东部明显的强降水(图9c2),预报范围和强度在第3—4小时明显增大、增强,至11—12时已接近实况(图9c4)。

图 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)
5.3 ZTD 3DVar对积分初期温度模拟效果的影响

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

图 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)
5.4 ZTD 3DVar对积分初期风场模拟效果的影响

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

图 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)
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资料同化有效增强了初始湿度场,对温、压、风初值无明显的调整,改进了积分初期水成物含量与分布位置进而改善了降水预报效果,修正了对照试验对辽宁东部地区强降水的明显漏报,且通过降水的反馈作用改进了温度与风场预报。

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

参考文献
陈敏, 范水勇, 仲跻芹等. 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 Meteor-ological 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. (in Chinese)
倪允琪, 周秀骥. 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. DOI:10.3321/j.issn:1001-8166.2006.09.001 (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 ground-based GPS ZTD in North China. Acta Meteor Sinica, 75(1): 147-164. (in Chinese)
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. DOI:10.1175/JAMC-D-15-0137.1
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. DOI:10.1175/MWR-D-11-00156.1
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. DOI:10.5194/angeo-27-2739-2009
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. DOI:10.1256/qj.03.26
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. DOI:10.1256/qj.04.102
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. DOI:10.1002/qj.693
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. DOI:10.1002/qj.1943
Dee D P. 1995. On-line estimation of error covariance parameters for atmospheric data assimilation. Mon Wea Rev, 123(4): 1128-1145. DOI:10.1175/1520-0493(1995)123<1128:OLEOEC>2.0.CO;2
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. DOI:10.1175/1520-0493(1999)127<1822:MLEOFA>2.0.CO;2
Desroziers G, Ivanov S. 2001. Diagnosis and adaptive tuning of observation-error parameters in a variational assimilation. Quart J Roy Meteor Soc, 127(574): 1433-1452. DOI:10.1002/qj.49712757417
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. DOI:10.1256/qj.05.108
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. DOI:10.3390/rs12030383
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. DOI:10.3402/tellusa.v38i2.11707
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. DOI:10.2151/jmsj.2004.507
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. DOI:10.5194/acp-17-13983-2017
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. DOI:10.1175/2007MWR2263.1
Mascitelli A, Federico S, Fortunato M, et al. 2019. Data assimilation of GPS-ZTD 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. DOI:10.1002/qj.2650
Mesinger F. 2008. Bias adjusted precipitation threat scores. Adv Geosci, 16: 137-142. DOI:10.5194/adgeo-16-137-2008
Parrish D F, Derber J C. 1992. The National Meteorological Center's spectral statistical-interpolation analysis system. Mon Wea Rev, 120(8): 1747-1763. DOI:10.1175/1520-0493(1992)120<1747:TNMCSS>2.0.CO;2
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. DOI:10.5194/amt-12-345-2019
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. DOI:10.1002/qj.2546
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. DOI:10.1016/S1464-1895(01)00091-6
Vedel H, Huang X Y. 2004. Impact of ground based GPS data on numerical weather prediction. J Meteor Soc Jpn, 82(1B): 459-472. DOI:10.2151/jmsj.2004.459
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. DOI:10.1002/qj.2661
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. DOI:10.1002/qj.508
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. DOI:10.1175/MWR-D-18-0418.1