1 引 言

湿度的时空分布变化很大，一直是数值预报初值场描述不够准确的要素之一，同时也是高分辨率短时降水预报不确定性的主要来源之一。地基全球定位系统(GPS)遥感技术根据大气折射导致的GPS卫星信号延迟量可提取大气水汽信息，为大气水汽监测提供重要的补充。从20世纪末开始，数值预报研究者就开始对地基GPS观测开展同化应用，以期改善分析场质量及预报技巧。

早期地基GPS观测在数值预报中的应用主要集中在其反演产品可降水量(PW)的同化应用上。国际上从20世纪末就开始了可降水量的同化应用研究(Kuo et al，1993，1996; Smith et al，2000; Guo et al，2005)。随着中国地基GPS气象观测网(GPS/MET)的建设进展，中国的可降水量资料同化研究也逐渐得到了重视(袁招洪等，2004；袁招洪，2005；张朝林等，2005；丁金才等，2007；楚艳丽等，2007，李昊睿等，2014，贝纯纯等，2016)。研究表明，同化可降水量可有效地改善大气低层湿度的初值质量，对提高短期数值预报的强降水预报性能有作用。随着可降水量同化应用对数值预报性能提高，尤其是降水预报技巧提高不断地被证明，中国业务预报中心开始了可降水量在数值预报业务系统中的同化应用(陈敏等，2010)。

由地基GPS测量大气水汽原理可知，GPS卫星发射带有时间标识的无线电信号，经大气层到达地面接收机时，由于大气折射的影响而产生延迟。地面接收机天顶总延迟(ZTD)包括流体静力学延迟和天顶湿延迟。静力学延迟是ZTD中的相对大量，与气压有很强的关系。天顶湿延迟则为ZTD中的相对小量。表征由于大气中水汽的存在造成的信号延迟量，需要通过基于包含地面气压和大气层平均温度等要素所建立的模型处理后才能得到对应的可降水量(Bevis et al，1994; Ware et al，1997)。这个处理模型是建立在要素观测和历史资料统计基础之上，而经验估算可能为所得到的可降水量导入额外的误差。比较而言，ZTD的误差来源比可降水量少，并由于更少涉及其他观测要素而更易于确定。逐渐地，地基GPS观测资料在数值预报中的同化应用研究更多集中于ZTD。相关个例研究(Vedel et al，2001; De Pondeca et al，2001; Cucurull et al，2004；Shoji et al，2011; Tilev-Tanriover et al，2014)表明，同化ZTD可比同化可降水量得到更好的效果。基于多个例或连续时段ZTD资料的同化试验(Vedel et al，2004; Poli et al，2007; Yan et al，2009；蔡雅婷等，2012)也表明，同化ZTD有助于改进对流层水汽初始场质量和降水预报性能。英国气象局自2007年就已经将ZTD业务同化应用到NAE(North Atlantic and European)数值预报模式中。

用于同化应用的观测资料往往是不完美的，可能含有由机器故障、记录出错或者资料处理和传输过程中的问题所造成的很大误差。含有这类误差的观测资料一般在数据中心观测资料收集整理过程中(面向观测资料的质量控制)会得到辨别和标识。此外，观测资料可能还包含有仪器误差、代表性误差等。诸如ZTD这种非直接观测，还会存在解算过程中由于算法精度导入的误差，这些误差可能是系统性的，也可能是随机的。

三维变分同化系统通过假设背景场误差和观测误差均服从高斯概率分布，由观测值、背景场以及它们的误差协方差矩阵来确定大气状态的最优线性无偏估计。即在观测资料进行同化应用前，除了开展面向观测资料的质量控制标注出“坏”数据，还需确定“好”数据的观测误差是否服从同化系统假设的观测误差无偏和高斯分布。如果不是，则需开展进一步的质量控制(面向同化应用的质量控制)，使质量控制后的观测资料尽量满足这一假设。由于真值无法获得，不能直接验证观测误差是否无偏和满足高斯分布，常用的可行验证方式是判断观测值与相应的背景场值的差值(观测背景差，OMB)是否满足高斯分布(Zou et al，2006; Qin et al，2010)。通过检查观测背景差大样本集的统计性质，诸如一阶矩均值，二阶矩标准差，三阶矩偏度系数以及四阶矩峰度系数，是否满足高斯分布所要求的特性，来确定观测背景差样本所对应的观测样本集是否满足同化系统关于观测性质的假设。通过面向同化应用的质量控制，使观测资料大体满足同化算法关于观测性质的假设，以助同化系统更为正确合理地同化分析观测数据所包含的大气信息，从而有效改进初值场质量。

如上述，欧洲自2000年前后就已经开始了ZTD资料的同化应用，也建立了相应的质量控制流程。比较普遍的做法(Poli et al，2007，Yan et al，2009，Bennitt et al，2012)是对逐站ZTD观测资料的观测背景差进行正态分布的假设检验，包括卡方检验(Chi-square test)，K-S检验(Kolmogorov-Smirnov test)等，剔除不能通过95%信度检验的站点；剔除出现异常值频率较高的站点，异常值通常通过不同处理中心解算值之间的差异来确定；剔除有效样本低于一定比例的站点，这一比例通常是40%—50%；剔除观测站点地形高度与模式地形差超过一定阈值的站点；对一个站点有多个处理中心解算数据的，选择观测背景差标准差最小的或者其频率分布最接近高斯分布特征的那个解算中心提供的数据。对通过上述步骤的各站点，将一段时间观测背景差均值作为偏差值，对观测进行偏差订正，时段通常为10 d以上。英国气象局(Bennitt et al，2012)的资料监视流程中确定不同处理中心提供的同一站点观测质量的标准是观测背景差标准差，标准差小于5 mm表征解算质量高，标准差大于10 mm 表征解算质量低。所有观测背景差大于5倍最大标准差的观测将被剔除，其研究表明观测背景差最大标准差的典型值为11 mm，即某ZTD观测资料的观测背景差大于55 mm，该观测即被剔除。这些质量控制策略在相应的研究中均表现出了其有效性。由于华北地区ZTD资料的解算软件和策略，有效数据率、异常值标准、解算质量标准，乃至资料同化观测算子都与上述研究大有不同，需要根据本地区ZTD数据特性来发展针对华北地区观测资料的质量控制策略和流程。

地处华北的北京、天津、河北和山西等省、直辖市已建立了较为完整的地基GPS观测网，能够提供有效观测的测站为146个。自2012年底开始，中国气象局北京城市气象研究所实现了对上述地区地基GPS观测资料的实时解算，数据时间分辨率为30 min。目前该观测资料还未得到有效应用。文中主要针对这些资料并基于北京市气象局数值预报业务系统BJ-RUC(BeiJing-Rapid Updated Cycling)，开展资料质量控制和同化应用研究，解决ZTD资料在数值预报系统中同化应用面临的问题。

2 天顶总延迟观测算子

本研究基于WRFDA(WRF model data assimilation system)三维变分同化系统对ZTD开展同化试验。在变分同化系统中，资料同化实际上是一个代价函数最小化的过程

(1)

式中，x为模式状态向量，x_b为从误差协方差为B的短期预报中获取的背景场变量，y_o为观测向量，H为观测算子，R为观测本身以及正演算子误差协方差(Lorenc et al，2000)。

ZTD的表达式为

(2)

式中，z为拔海高度，N为大气折射率(以10⁶的比例因子被放大)(Smith et al，1953)，可以表达为

(3)

式中，p_d、p_v分别为干空气气压和水汽压，T为温度，k₁、k₂、k₃>分别为空气折射率常数(Bevis et al，1994)。将状态方程p_d=ρ_dR_dT,p_v=ρ_vR_vT分别代入式(3)，可得

(4)

将p_d=p-p_v，以及水汽压定义p_v=代入式(4)，可得

(5)

式中，q为比湿，结合虚温定义T_v=T(1+)，式(5)可表达成

(6)

在模式中对大气折射率进行积分时，采用了大气折射率随高度指数减小的假定

(7)

从而避免了假定两模式层间大气折射率为常数造成的对ZTD低估的问题(Bennitt et al，2012；Poli et al，2007)。模式格点上背景场ZTD表达式为

(8)

式中，

(9)

(10)

式(9)为模式格点上从模式顶(k=k_top)至模式底(k=1)积分得到的背景场ZTD，其中，k为模式层；式(10)为模式格点上模式层顶至卫星的ZTD估计(Saastamoinen，1972)，p_top为模式顶层气压，z_top为模式顶层高度，φ为格点所在位置的纬度。使用背景场中测站位置周围4个格点上ZTD进行线性插值，即可得到该测站位置的背景场ZTD。

由式(2)可知，ZTD是从卫星高度至接收机高度(通常取测站高度)的积分量。观测算子分两步求解与观测ZTD相匹配的背景场ZTD。第1步采用上述各式求解测站位置从卫星高度到模式地形高度积分得到的天顶总延迟，记为ZTD_m。由于数值模式中地形分辨率有限，测站往往并不恰好位于模式地形上，尤其在地形复杂区域。测站高度可能低于模式地形，也可能高于某个模式层高度，第1步计算得到的ZTD_m与测站观测ZTD并不匹配。所以，第2步需要计算由于模式地形和测站高度不一致导致的背景场ZTD订正值，通常称为地形差补偿值d_ZTD。

对于测站高度低于模式地形高度的测站，地形补偿值为以测站高度为下限、模式地形高度为上限计算得到的背景场ZTD，该测站位置上的地形差补偿值为

(11)

式中，h_m为模式地形高度，h_o为测站高度，N₁为模式最低层上的模式大气折射率。

对于测站高度高于模式地形高度的测站，地形补偿值则为以模式地形高度为下限、以测站高度为上限积分得到的背景场ZTD。该测站位置上的地形差补偿值为

(12)

式中，h_n为高于测站高度且最接近测站高度的模式层高度，N_n为该高度上的模式大气折射率，h_o为测站高度，N_z为低于测站高度的各模式层上模式大气折射率。

订正后的背景场ZTD表达式为

(13)

式中，ZTD_m可从卫星高度到模式地形高度积分而得到，地形订正非常重要，尤其是在站点高度与模式地形差别较为明显的情形下，表 1分别给出了站点高度高于和低于模式地形的站点相关信息进行定量说明。以北京斋堂站为例，该站的站点高度为780 m，3 km分辨率的模式地形为444 m，地形差高达336 m。2014年5月8日08时(北京时，下同)，该站ZTD观测为227.9 cm，由式(1)—(8)计算出美国环境预报中心(NCEP)GFS(Global Forecast System)分析场中相应位置上的ZTD_m为216.9 cm，观测值与ZTD_m之差为11 cm。由式(9)—(10)计算GFS分析场中相应位置上由于模式地形低于站点高度造成的补偿值d_ZTD为-9.6 cm，因此，最终GFS分析场中斋堂站点上观测背景差为1.4 cm。同化系统通过观测背景差来引入观测带来的大气信息，尽管地形补偿值-9.6 cm相对于观测值227.9 cm是一个小值，但相对于观测背景差11.0 cm却是一个量级相当不可忽略的值，可见地形订正在ZTD观测算子中的不可或缺。

文中主要针对华北地区的地基GPS观测展开，目前这一地区能够有效提供观测数据的地基GPS测站为146个，大致均匀分布在京津冀晋4省市。其中北京站点密度较高，河北北部站点密度略低(图 1)。华北地区位于中国地势第二阶梯和第三阶梯交界区域，自西北至东南从山脉到平原，测站高度从高于海平面1992 m至低于海平面3 m，跨度较大。由于模式水平网格分辨率的限制，模式地形与测站高度差异明显。在水平分辨率3 km的模式网格中，35%的测站位置上模式地形与测站高度的差超过100 m，3%测站位置上模式地形与测站高度的差超过300 m。由此可见ZTD观测算子中地形补偿的必要性。如果所使用的同化观测算子中没有有效的地形订正算法，需要将模式观测地形差达到一定阈值的站点加以剔除，以免由于过大地形差引入误差影响同化效果。使用不同同化系统的不同研究人员对这一阈值有不同定义，Yan等(2009)定义该阈值为150 m，而Bennitt等(2012)则定义该阈值为300 m。

3 ZTD观测及其质量控制 3.1 ZTD观测特征

ZTD包含大气水汽信息。单站ZTD观测的变化主要由该站天顶方向大气柱中水汽含量和地面气压的变化决定。图 2a给出了北京琉璃河站2014年ZTD和可降水量观测时间序列，可以看到ZTD的变化基本上反映了可降水量，即测站天顶方向大气柱中水汽含量的变化。7—8月(夏季)华北地区水汽丰沛时期ZTD高达270 cm且变化幅度很大，而12—1月(冬季)华北地区干燥时期则低至235 cm左右，且变化幅度较小，春、秋两季无论是ZTD还是其变化幅度均基本居于冬、夏两季之间。图 2b给出了2014年北京琉璃河站ZTD观测背景差时间序列，其中背景场选用GFS分析场。ZTD观测背景差总体为-2—2 cm，同时表现出低频波动和尖锐凸出。低频波动基本反映出随着天气系统的生消、发展和移动，观测与背景场差异的变化趋势，且低频波动的振幅表现出季节特征。尖锐凸出表明某些时刻观测与背景场之间差异巨大，可能是观测资料包含较大误差造成的，也可能是观测资料反映真实大气状况而背景场受限于空间分辨率不能很好地描述水汽的高频变化造成的。这种情况下观测背景差所对应的ZTD观测往往不会被同化系统吸收，即便其表达的是真实大气状况，因为与背景场偏离过大的观测信息的吸收可能会造成背景场中物理量空间不连续或者物理量间不平衡进而破坏后续积分过程。通常同化系统会包含简单的质量控制方案，即观测误差倍数方法(观测背景差绝对值超过观测误差的一定倍数，通常取5倍或者3倍，则剔除观测背景差所对应的观测)。但观测误差，尤其是一些局地新类型观测资料的观测误差的估计也还存在一定问题，基于观测误差开展质量控制这一方式有其自身的不足之处。因此，研发合理和适用的ZTD观测质量控制方案是该数据在科研和业务中进行同化预报应用之前不可忽视的重要环节。

3.2 观测质量控制

研究收集了2013年1月—2014年12月逐日02、08、14、20时的北京、天津、河北和山西四地ZTD观测数据，建立一套ZTD观测序列。选用GFS分析场作为背景场，基于WRFDA同化系统，建立了ZTD观测序列对应的ZTD观测背景差序列。质量控制方案的研发和合理性验证均基于上述样本序列开展。

3.2.1 质量控制方案

对样本序列开展统计性质检验的一个基本前提是序列中样本满足“独立同分布”的假设，即样本之间是相互独立的，序列内各样本服从同一个分布。由图 2a所示的由2014年ZTD观测序列分布可知，ZTD在不同季节、甚至不同月的分布特征都不尽相同，离散程度有显著差异。当序列由单月内ZTD观测样本组成，则该序列中的样本更接近“独立同分布”的假设。因此，文中基于单月样本组成的序列开展质量控制，分析质量控制前后序列的统计特性，来确认质量控制方法的适用性，并定义不同序列为：①ZTD观测的单站单月序列SSobs(Single station Single month obs serial)，序列中样本为SSobs_i(i=1，2，…，n)，ZTD观测背景差的单站单月序列SSomb，序列中样本为SSomb_i(i=1，2，…，n)，n为单月样本数(当月天数×4个样本/d)。②ZTD观测的所有站单月序列ASobs(All station Single month obs serial)，序列中样本为ASobs_ij(i=1，2，…，n;j=1，2，…，m)，ZTD观测背景差的所有站单月序列ASomb，序列中样本为ASomb_ij(i=1，2，…，n;j=1，2，…，m)，m为站点数。

不同质量控制环节基于上面定义的不同序列，针对不同的观测性质进行质控，将不能通过质控的观测标记为可疑观测，仅保留通过样本进入下一个质量控制环节。

(1) ZTD的解算误差和稳定性检查(QC1)

GPS卫星信号传输延迟量的确定取决于诸多因素，解算过程将这些因素的不确定性综合表达为延迟信号的解算误差，即解算软件使用模型和算法估计ZTD过程中的估计误差，该值越大表示解算值越离散，可信度越低。中国气象局北京城市气象研究所解算的华北地区ZTD观测大样本序列的统计分析表明，当ZTD观测对应的解算均方根误差大于1.0 cm时，该数据通常不够可靠。因此，QC1将解算均方根误差大于1.0 cm的ZTD观测标记为可疑。此外，一方面由于文中的质量控制方案基于序列的统计性质开展，当序列中样本数过少时，统计量可能不足以表达序列性质；另一方面测站若缺测频率较高也一定程度上表明该测站不稳定。因此，QC1将SSobs序列中实有样本数少于应有样本数1/3的站点标记为可疑站点。以1月为例，某站点的SSobs序列应有样本数为124个(31×4)，如果该站该月序列少于37个，则该站点标记为可疑站点。由QC1标记出的可疑观测和可疑站点的观测，均不参加后续环节的质量控制。

图 3为2014年5月华北地区各站点SSobs序列进行QC1的情况。可以看到，标示为红色的测站为有效样本少于应有样本数的1/3的测站。这些测站的ZTD观测不能进入后续质控环节，进而无法参加同化应用。标示为蓝色的测站为通过QC1的测站，这些测站上解算均方根误差小于1.0 cm的ZTD观测将进入后续质控环节。图 4展示了2013—2014年逐月ASomb序列在QC1前后偏度系数的变化。由图 4可知，剔除了解算误差过大的观测以及不稳定站点的观测，由ZTD观测背景差组成的ASomb序列的偏度系数在一些月份中有明显变化，更为接近高斯分布偏度系数0，尤其以2013年3月为明显。表明剔除不稳定测站和解算误差过大的观测确实有助于序列更加接近同化系统对资料性质的假设。

(2)观测背景差的双权重检查(QC2)

传统方法在确定样本序列中的离群值时，通常基于序列中所有样本(A_i，i=1，2，…，n)估计序列均值m和标准差σ，将偏离均值达到标准差的某个倍数的样本确定为离群值，通常取3倍，即满足的样本为离群值。但是序列中离群值的存在，往往会使得传统方法计算得到的均值和标准偏差被一定程度地歪曲，进而影响离群值的有效确定。Lanzante(1996)的研究表明使用双权重检查方法估计均值和标准差，较之使用传统方法估计均值和标准差更少受到离群值的影响。

在双权重估计方法中，对于由n个样本组成的序列(X_i，i=1，2，…，n)，中值为M，序列中每个样本与中值M的差值绝对值的均值为X_MAD。每个样本(X_i)的权重函数定义为

(14)

式中，c为一个阈值，所有偏离中值达到一定距离的数据均给定0权重。文中这个阈值取7.5。此外，在权重函数的计算中|w_i|>1时，则设置|w_i|=1。

双权重均值X_bi和双权重标准差σ_BSD分别定义为

(15)

(16)

从式(15)、(16)的定义可以看到，双权重均值X_bi和双权重标准差σ_BSD的计算中，位于序列分布中心的数据要比位于序列分布边缘的数据得到更大的权重。双权重方法计算出的均值和标准差对于序列分布边缘数据更不敏感，相对于传统方法给处于序列分布不同位置的样本相同权重的做法，双权重方法在计算均值和标准差时更为合理。

X_bi和σ_BSD用于计算每个样本的Z评分

(17)

Z评分用做质量控制标记的判据。当一个样本的Z评分大于设定阈值时，该样本则被标记为所在序列中的离群值。文中Z评分阈值设定为3，即针对SSomb序列进行QC2，将序列中Z评分大于3的样本确定为离群值，被确定为离群值的ZTD观测背景差样本所对应的ZTD观测将被标记为可疑观测，可疑观测不进入后续质量控制以及同化应用。

从2014年5月北京琉璃河站的SSomb序列频率分布(图 5)可知，该序列中多数样本集中分布在-2—3 cm，少数小于-4 cm或大于6 cm的样本偏离样本主体位于分布函数的尾端。分别使用传统方法和双权重方法确定的该序列非离群值判断条件和范围如表 2所示。可以看到，两种方法确定的序列均值和标准差明显不同，尤其是标准差。受位于序列分布边缘数据的影响，传统方法计算得到的标准差偏大，而双权重方法由于给予序列分布边缘数据更小的权重，因而受其影响较小，计算得到的标准差偏小。因而两个方法确定的非离群值范围也显著不同。由图 5可以看到，传统方法表明该序列中无离群值，而双权重方法则确定出序列中大于4.83 cm的4个样本和小于-2.98 cm的3个样本为离群值。实际上基于序列频率分布，从视觉直观判断也能得出与双权重方法类似的结论，但传统方法却受方法本身的限制而无法做出正确判断。

图 6为2014年5月的ASomb序列的散点分布，其中红色点为QC2确定的SSomb序列的离群样本，大多偏离ASomb序列的对角线位置，其所对应的ZTD观测被标记为可疑观测，不进入后续质控和同化应用。由于针对SSomb序列进行QC2只能甄别出单站序列中的离群样本，而相当一部分测站的SSomb序列本身有或正或负、或大或小的系统偏差，因而尽管图 6中黑色点由通过QC2的SSomb样本组成，仍然有样本处于明显偏离对角线的位置。这说明还需要有其他质量控制手段来进行进一步的检查和控制，例如针对SSomb序列进行系统偏差订正。

(3) 偏差订正(QC3)

图 7为华北地区各站点2013—2014年ZTD观测背景差序列的均值和标准差。很明显相当多测站的值存在系统偏差。QC3将对各站点SSomb序列进行偏差订正，偏差订正值取各站点SSomb序列的均值。

(18)

式中，i=1，2，…，n，为站点序号，j=1，2，…，m，为单站的样本序号，O_Meani为第i个站点的SSomb序列均值，O_i，j为第i个站点偏差订正前的SSomb序列第j个样本，O_correcti，j为第i个站点偏差订正后的SSomb序列第j个样本。

针对SSomb序列开展偏差订正，一方面是由于不同测站系统偏差情况各异，需要针对不同情况个别处理。另一方面取月序列均值能表达观测和相应背景场的系统偏差，偏差的高频噪声随着时间平均会相对平滑，同时具有较好的代表性。QC2已经标记出SSomb序列的离群值，避免了离群值存在造成的系统偏差估计的错误。QC2与QC3共同完成了对SSomb序列系统偏差的合理估计和相应订正。图 8为QC3前后ASomb序列的偏度系数。很明显，QC3使得ASomb序列的偏度系数明显更为接近高斯分布偏度系数值0。

(4) 观测背景差的标准差检查(QC4)

QC1—QC3各个环节都是基于SSobs或者SSomb进行不同性质的检查，QC3完成后，有些站点的SSobs和SSomb序列样本数少于应有样本数的一半，序列的相应统计特征明显区别于其他站点，尤其是SSomb序列的标准差。分析ASomb序列统计特性发现，SSomb序列标准差明显偏大的站点，往往致使其所在的ASomb序列峰度系数异常偏大。基于这个发现，确立了针对ASomb序列的QC4环节：分别计算各站点SSomb序列和同月ASomb序列的标准差，当某站点SSomb序列的标准差大于同月ASomb序列标准差且达到既定阈值，则该站点被标记为可疑站点，可疑站点的ZTD观测不参加同化应用。文中这个阈值确定为ASomb序列标准差的1.25倍。

图 9为2014年5月通过QC3的华北地区所有站点的SSomb序列的标准差，其中绿色直线对应的值为ASomb序列的标准差。可以看出，不同站点的SSomb序列标准差有明显差异，部分站点的SSomb序列标准差与ASomb序列标准差差异很大。红色柱表示该站点的SSomb序列标准差超过ASomb序列标准差的1.25倍，未能通过QC4，这些测站的ZTD观测将被标记为可疑观测，不参加同化应用。由QC4前后ASomb序列的峰度系数对比(图 10)可以看出，将QC4标记出的可疑观测排除后，ASomb序列的峰度系数明显更为接近高斯分布的峰度系数值3。

3.2.2 质量控制方案的合理性

QC1针对SSobs序列，QC2—QC3针对SSomb序列，对序列中可疑样本进行标记和剔除，对序列的系统偏差进行订正。QC4针对SSomb与ASomb序列标准差比例，对可疑测站进行标记和剔除。通过这一套流程，对ZTD观测开展面向同化应用的质量控制，实现质控后的观测资料大体满足观测误差无偏且高斯分布这一变分同化系统的假设。通过质控前后序列频数分布和统计特征可以确认质量控制方案的合理性。图 11为2014年5月ASomb序列在质量控制前后的频数分布。质量控制前序列频数分布跨度很宽，在远离样本主体的区间内仍有样本分布，质控后这些离群值都被剔除，使序列分布跨度明显减小。质量控制前序列在均值附近分布频数显著高于拟合的高斯分布，在1—2倍标准差区间内分布频数则明显低于拟合的高斯分布。与之相比较，质量控制后序列除了在均值附近频数分布高于拟合高斯分布外，在整个区间上均很好地贴合了拟合高斯分布。表 3给出了该序列在质量控制前后的统计特征量，在标记并剔除了7.5%样本的情况下，使序列从质量控制前的负偏差、负偏度系数和过大峰度系数，变为质控后的0偏差和接近高斯分布要求的偏度系数和峰度系数，标准差也一定程度减小。基本达到了同化系统对观测资料的性质假设，表明文中所建立的质量控制方案的合理性。

图 12展示了针对2013—2014年华北地区ZTD观测资料，使用文中所建立的质量控制方案完成质控后，序列的统计特征和质控剔除数据比例。图 12a为质量控制后ASomb序列的均值和标准差，可以看到均值为0，满足无偏分布，标准差基本呈夏季高、冬季低、春秋居中的特征。图 12b为质量控制后ASomb序列的偏度系数和峰度系数，可以看到偏度系数基本在0左右而峰度系数基本在3左右，大体与高斯分布统计量特征一致。与此同时，质量控制实现这样的效果并未过多牺牲观测数据的使用效率，由图 12c所示的质控中被标记和剔除样本数的百分比可以看到，除个别月份标记和剔除样本数占质控前总样本的14%左右，其余月份剔除率均为2%—8%。

此外，文中还使用探空观测诊断ZTD，对质控前后的ZTD观测进行了验证。选取北京探空站(站号54511)，以及距离54511站最近且拔海高度相近的两个GPS站点：大兴气象局站(DAXN)和朝阳气象局站(CHAO)。站点信息分别为：北京探空站(39.80°N，116.47°E，33 m)，大兴气象局站(39.72°N，116.35°E，40 m)，朝阳气象局站(39.95°N，116.49°E，30 m)。选取54511站2014年6月1日—8月31日，每天3次探空(08、14和20时)观测诊断ZTD为真值，计算该时段内GPS站点观测ZTD的平均偏差和均方根误差(表 4)。可以看出，质量控制后两个GPS站点观测ZTD误差均值明显减小，均方根误差也明显减小，表明质量控制方案中偏差订正是有效的。

3.2.3 业务应用中质量控制策略

以上质量控制方案各环节均基于单月序列开展，由于所针对的样本为历史观测数据，所有单月序列由一个自然月内所有观测样本组成。在面向业务数值预报同化应用的实时质控中，这一序列组建方式显然不可行，因此，需要在实时应用中调整策略。一种可行的方式是：以当前时间为截止时间，以距离当前时间30 d的对应时间为起始时间，取这一时段内所有观测数据组成单月序列，并基于这一序列开展质量控制。当前时间观测若被质控标记为可疑观测，则不参加当前时次业务同化应用。

4 ZTD观测的数值预报同化应用 4.1 数值预报模式和试验设计

基于北京市气象局数值预报业务系统BJ-RUCv2.0开展地基ZTD同化和预报试验，BJ-RUCv2.0系统预报区域为两层嵌套，分辨率分别为9和3 km。根据背景场来源分为冷、热启动预报两类，其中冷启动预报为第1个预报循环，该循环以GFS预报场作为同化背景场，而热启动预报则在上次循环的预报场基础上进行同化分析。循环系统由逐日08时冷启动开始，其后每隔3 ｈ进行一次三维变分更新同化，即在当日11、14、17、20、23时和次日02、05时，同化当时次观测资料并进行预报。由于GPS观测集中在北京、天津、河北和山西4个地区，同化预报试验均在BJ-RUCv2.0的3 km分辨率区域(图 1)开展，GPS观测覆盖区域占整个预报区域1/4左右。

数值试验设计如表 5所示。数值试验时段为2014年5月31日—6月11日共11 d，每天8次预报共88次预报试验，每个循环的预报时效均为24 h。这11 d内，北京及周边地区出现了7次明显的对流性降水，大多伴有短时大风和雷电。通过比较各试验的降水预报性能，确认文中所建立的ZTD质量控制方案的有效性，确认同化ZTD对降水预报性能的影响。在同化ZTD试验中，ZTD观测误差取为1.25 cm，观测误差的调整将在今后进一步开展评估。

4.2 ZTD观测同化应用对降水预报的影响 4.2.1 ZTD同化对降水预报的整体性能影响

评估针对整个预报区域试验时段中4个试验开展的逐3 h降水预报，所用实况资料为预报区域中所有中国国家级常规地面站的3 h降水观测。模式网格点降水预报通过双线性插值到地面观测站点位置生成模式的测站预报值。评估采用国际上通用的TS(Threat Score)评分(Donaldson et al，1975)。

图 13给出各试验的5 mm/(3 h)、10 mm/(3 h)两个降水阈值的TS评分。总体上，可降水量、ZTD_NOQC和ZTD_QCED试验的预报性能优于CNTL。比较可知，除5 mm阈值的0—3 h预报时效外，ZTD_QCED试验TS评分在其他预报时效上均明显优于ZTD_NOQC试验，证明了该质量控制方案的合理性；ZTD_QCED试验TS评分在所有预报时效上优于或相当于可降水量试验，表明ZTD观测同化应用对于降水预报性能的提高作用比可降水量观测同化应用要更为明显，为ZTD资料在业务预报同化应用中取代可降水量资料提供了科学依据；除5 mm阈值的0—3 h预报时效和10 mm阈值3—6 h预报时效外，ZTD_QCED试验TS评分在其他各预报时效上均优于CNTL试验，证明了同化质量控制后的ZTD观测能够有效地提升业务系统的降水预报性能。

20 mm/(3 h)、30 mm/(3 h)阈值降水的TS评分(图略)均在0.1以下，甚至0.01以下，各试验也无明显优劣之分。40 mm/(3 h)及以上阈值降水在试验时段内仅在2个3 h累计时段中出现过，也基本无预报技巧。由于GPS资料的同化能影响分析时次水汽分布并进而影响不稳定条件，但如果控制试验对降水的预报能力本身很弱，同化GPS，不论可降水量还是ZTD，以及是否质量控制都很难产生明显的正面影响。

TS评分从出现降水事件的站点数中预报正确的站点数所占比例的角度来描述不同试验的降水性能，下面尝试直接比较不同试验的10 mm/(3 h)阈值降水在各预报时效上的预报准确站数和预报错误站数，来细致直观认识ZTD资料参加同化对降水预报的影响。

众所周知，探空观测是可靠和有效的高空大气信息来源，但其时间分辨率有限。文中开展的快速更新同化预报数值试验中，有探空观测参加同化的只有08、14、20时启动的预报，11、17、23、02、05时启动的预报则无探空观测可同化应用。图 14a、b分别给出了有探空观测参加同化的预报平均的预报正确站数和预报错误站数。图 14a中3个试验在不同时效上的预报正确站数各有不同，总体上并无哪个试验有明显优势。图 14b中ZTD_NOQC和ZTD_QCED试验预报错误站数在多数预报时效上都小于CNTL试验，可见同化ZTD资料可以有效地减少预报错误站数，即有效改善空报、漏报现象。ZTD_QCED试验在所有预报时效上预报错误站数都最少，表明质量控制后的ZTD资料同化使得这种改善效果更为明显。图 14c和d为无探空观测参加同化的预报平均的预报正确站数和预报错误站数。图 14c中大多数时效上ZTD_QCED试验的预报正确站数最多，图 14d中大多数时效上ZTD_NOQC和ZTD_QCED试验的预报错误站点数少于CNTL试验，且预报错误站数比CNTL试验偏小的程度比图 14b中更显著。表明在无探空观测参加同化的预报中，同化ZTD资料改善空报、漏报的作用更为显著。同化质量控制后的ZTD在更为显著地改善空报、漏报现象的同时，还有效地增加了预报正确站数，对预报性能的提升作用更为凸出。可见ZTD资料在无探空观测参加同化的预报时次中，作为丰富高空大气信息的有效来源，对改善初值场质量有明显作用。也从资料同化对预报性能的影响角度证明了本研究所建立的质量控制方案的合理和有效。

由于目前试验中同化GPS观测站点仅限京津冀晋4省市，预报区域中相当大区域并无GPS观测站点，用以改善预报区域湿度初始条件的观测信息有限。预期覆盖更大区域的GPS观测用于同化能进一步提升降水预报性能。

4.2.2 ZTD在个例预报中的同化应用分析

2 014年6月6日午间至傍晚，京津冀地区出现了一次短时降水过程，强对流云团自北向南影响京津冀大部分地区，其中河北出现了雨强达58 mm/(3 h)、北京、天津均出现了雨强大于35 mm/(3 h)的短时降水。图 15a为6月6日14—17时3 h累计时段内降水实况，达到10 mm/(3 h)及以上雨强的降水区域主要位于河北东北部、天津地区，成带状分布，京津冀其他地区大多为小雨。图 15b、c为6月6日14时循环启动的CNTL和ZTD_QCED试验在14—17时的3 h累计降水预报。CNTL试验基本上模拟出了雨带的走向和强度，但北京东南地区20 mm/(3 h)强度降水出现了较为明显的空报，与此同时对天津南部20 mm/(3 h)强度降水出现了较为明显的漏报。此外，河北北部出现20 mm/(3 h)强度的降水区域与实况相比也偏大，河北南部10 mm/(3 h)强度降水出现空报。ZTD_QCED试验则较好地预报出天津南部的20 mm/(3 h)强度降水，有效地减弱了CNTL试验中在北京东南部的空报和天津南部的漏报现象。ZTD_QCED试验还有效地减弱了河北南部10 mm/(3 h)降水的空报现象。

在与北京东北部和天津北部相邻的河北地区(40.5°N，117°E)，ZTD_QCED试验与ZTD_NOQC试验相比出现了20 mm以上量级的空报。究其原因，质量控制方案是建立在统计性质分析基础上的，从统计角度的最优未必是具体某个观测的最优，因此质控序列中某个观测可能处在序列边缘，对预报产生了负面影响。但总体来看，质量控制后的观测同化应用后得到了更优的预报性能，显示目前质控方案的合理性，同时也表明还需要开展更深入分析，进一步优化质量控制方案。

图 16为ZTD_QC与CNTL试验分析场模式第7层(约900 hPa)比湿差值和最大对流有效位能差值分布。与CNTL试验相比，同化ZTD观测资料后，分析场中天津北部和南部有很明显的增湿现象；北京、天津间的河北廊坊出现了南部湿度增大，北部湿度减小的现象。最大比湿增加(减少)量达4 g/kg(2 g/kg)。在湿度增大(减小)的位置，最大对流有效位能也对应地表现为增大(减小)。比湿和最大对流有效位能增大(减小)的中心位置与图 15b中CNTL试验漏报(空报)位置大致对应，表明ZTD资料的同化导致的模式大气湿度调整是合理有效的，有利于形成更接近真实大气中局地对流形成和发展的环境。

5 结论和讨论

针对2013—2014年北京、天津、河北和山西的地基GPS天顶总延迟观测，发展了面向数值预报同化应用的ZTD质量控制方案。使用3 h快速更新循环数值预报平台进行了有无同化GPS观测、同化不同GPS观测的平行试验。通过观测资料集在质量控制前后的统计特性对比分析和平行试验预报性能检验评估，对所发展的ZTD观测质量控制方案性能以及ZTD观测同化应用对预报性能的改善有了初步认识：

(1) 文中发展的ZTD资料质量控制方案从测站稳定性和解算精度、离群值、偏差订正以及异常标准差等考察点切入，从不同角度检查和标记出致使观测序列统计特性偏离高斯分布特征的样本。应用该方案后，ZTD观测序列统计特性更加接近三维变分系统对观测资料的统计性质假设，表明文中所建立的质量控制方案的合理性。

(2) 从2014年5月31日08时起开展了为期11 d 的连续3 h快速更新循环预报试验，5和10 mm阈值降水预报效果评估表明，ZTD_QCED试验TS评分多数时效上优于ZTD_NOQC试验，证明了本研究所建立的质量控制方案的有效性，ZTD_QCED试验评分多数时效上优于CNTL试验，表明ZTD同化应用可以有效地提升预报系统的降水预报性能；ZTD_QCED试验在所有预报时效上优于或相当于可降水量试验，表明ZTD观测代替可降水量在数值预报业务中开展同化应用的可行性。

(3) 对比3 h快速更新循环数值试验中有无探空观测参加同化的预报评估结果可以看到，同化ZTD资料能改善CNTL试验的空报、漏报，同化质量控制后的ZTD使得这一改善更为明显。在无探空观测参加同化的预报中，同化质量控制后的ZTD观测对预报性能的提升更为凸出。说明ZTD资料是对高空大气信息的重要补充，对改善初值场质量有明显作用。

(4) 2014年6月6日京津冀地区降水个例的同化预报分析表明，同化ZTD资料为初值场提供了更为细致的湿度调整信息，使得同化ZTD后的分析场中对流发生、发展条件更为合理。

总体来看，本研究发展的ZTD资料质量控制方案是合理有效的，使得质量控制后的资料统计特征更加满足变分系统对资料的性质假设。质量控制后的ZTD同化应用有效地提高了数值预报的局地降水预报性能，特别是在无探空观测的预报时次。文中同化预报试验中ZTD观测误差采用了主观设定值1.25 cm，下一步工作中需要使用合理方法对这一主观设定值进行优化调整或者采用更为科学的方法另行估计。