水成物对GNOS掩星弯曲角同化的影响评估

袁炳， 希爽,2， 廖蜜,2， 马刚,2*， 张鹏,2

1 国家卫星气象中心， 北京 100081 2 中国遥感卫星辐射测量和定标重点开放实验室， 北京 100081

0 引言

全球导航卫星系统(GNSS，Global Navigation Satellite System)掩星探测是20世纪80年代后期开始兴起的探测地球大气的新方法.通过测量穿过地球大气层时，由于地球大气层的温度、湿度和压力所引起的GPS信号延迟(由减速和弯曲引起)，来获得地球大气的温度、湿度和压力信息.从全球分析来看，掩星资料的反演精度与常规探空数据相当，甚至可作为探空数据的替代产品(Reigber et al.，2003).掩星资料可有效应用于气候分析(罗佳等，2018)，甚至能有效应用于高层云参数的反演(严卫等，2012)，而其更为广泛的应用则是为数值天气预报同化系统提供新型的大量观测数据.

早期掩星资料的同化通常采用Abel变换将弯曲角转换为射线切点上的折射率，观测算子在模式空间把变量内插至切点位置，给出在模式空间上切点处的折射率.研究表明，折射率资料同化对于重构水汽垂直廓线、改进温度场十分有效(Zou et al.，1995；Kuo et al.，1998；王鑫等，2005)，但也带来一定反演误差(宫晓艳等，2007).一般掩星连线与大气层的切线长度大于300 km，由于局地球面假设，Abel变换在对流层低层会导致系统性负偏差(Kursinski et al.，1997；Kursinski et al.，2000；Sokolovskiy，2003),此时在Abel变换中需利用沿切线方向的大气密度进行订正，但这仍不能补偿切线方向大气变量非均一性带来的正演误差.射线追踪非局地观测算子进一步改进了弯曲角资料同化效果(Zou et al.，1999，2002).利用GPS L1和L2双频信号组合进行电离层订正后的弯曲角是GPS掩星观测中较为原始的观测量.弯曲角资料的同化避免了折射率同化方案中引入的反演误差，且弯曲角误差特性相对简单.在目前的 ECMWF(European Centre for Medium-Range Weather Forecasts，欧洲中期天气预报中心)业务同化GPS掩星数据中，同化方案即为同化GPS掩星弯曲角资料.业务经验表明，弯曲角资料同化对对流层上层和平流层下层的温度初始场有明显改善，并在预报中对500 hPa距平相关的提升有重大贡献(Healy，2008；Poli et al.，2010；Rennie，2010).NCEP(National Centers for Environmental Prediction)的经验进一步指出，掩星弯曲角资料的同化对数值预报的贡献远大于折射率的同化(Cucurull et al.，2013).

在全球数值预报中，采用非局地观测算子进行弯曲角同化，理论上比采用局地观测算子有优势，而在区域数值模式中，因模拟区域的缩小、模式顶高度的下降和侧边界的影响，通常仅采用局地观测算子.然而，随着近年来计算条件的提高，区域模式的模拟区域增大，分辨率极大提高，中小尺度信息丰富，基于球对称假设的局地观测算子逐渐不满足理论需求.因此，非局地掩星弯曲角观测算子在区域模式中得到了一定应用.一些研究表明，无线电掩星弯曲角资料射线追踪算法同化对于改进台风个例的路径模拟误差较为助益(毕研盟等，2013).

云雨对卫星探测技术有不同程度的影响，大气中低频微波辐射传输受到降水和大粒子的云水散射影响，高频微波辐射传输受到云冰等冰相粒子的散射影响，中频微波辐射受到云水粒子的散射影响(Kuo et al.，2004；Ma et al.，2011；Bi et al.，2012；Legrand et al.，2016；Geer et al.，2017).而通常认为GPS掩星探测受云和降水影响较小，在观测算子正向模拟计算中可忽略不计.然而近年的研究表明，GPS反演大气廓线在云区存在系统性误差，对温度为正，对水汽为负(Zou et al.，2012)，当云水含量大于0.6 g·m-3时，水相和冰相云粒子将对GPS信号产生影响.同时，折射率在对流层低层存在负的系统性偏差，对于不同云型，正演折射率的影响在1.2%左右(Lin et al.，2010).尽管在全球模式同化掩星折射率时忽略了云水的影响(Rocken et al.，1997)，但全球平均云覆盖达60%左右，全球模式中掩星信号显示了对云水的弱敏感性(Schmidt et al.，2005).

FY-3c新增了全球导航卫星掩星探测仪(GNOS，Global Navigation Satellite System Occultation Sounder)，实现了北斗导航系统(BDS，BeiDou Navigation Satellite System)和全球定位系统(GPS)双系统兼容掩星探测，大大增加了掩星观测次数.平均每天在轨接收北斗掩星事件约200次，GPS掩星事件500多次(王树志等，2015)，为数值天气预报和全球气候分析提供了更多原始数据.国家卫星气象中心(NSMC, National Satellite Meteorological Center)的业务经验表明，GNOS反演大气温湿廓线的精度与国外GPS的反演精度相当(Yang and Zou，2012).

本文中基于区域三维变分同化系统WRFDA (WRF model data assimilation system)，构建了FY-3c GNOS弯曲角的同化算子，针对台风天气系统开展FY-3c GNOS掩星弯曲角资料的同化试验.其中，针对云型丰富、云水含量空间梯度变化复杂的特点，在同化系统中引入了云水、云冰等多相态水成物对GNOS正演弯曲角的影响，评估了区域数值模式中水成物对掩星观测同化和预报的影响.本文分为四部分，第一部分为引言，介绍了当前掩星资料同化和风云卫星掩星探测技术的最新进展，第二部分为理论和方法介绍，重点描述了射线追踪算法和水成物对掩星弯曲角正演的影响，第三部分为试验和分析，利用台风个例评估了GNOS掩星弯曲角同化对台风天气模拟的影响，以及水成物对掩星弯曲角资料同化的影响.第四部分为全文总结.

1 理论和方法

1.1 掩星弯曲角射线追踪算法

从已知的GPS卫星、LEO卫星的精密星历以及大气状态，利用射线打靶法可以获得信号的轨迹，从而得到弯曲角和影响参数序列，完成三维追踪过程.但其计算量较大，资料的实时处理存在困难.为此，Zou等(1999)发展了掩星平面射线追踪算子.Liu和Zou(2003)又改进了射线追踪算子的精度和效率.无线电掩星观测资料的水平分辨率在大气的中低层通常大于100 km且分布不均匀，而射线追踪法是当前国内外解决无线电掩星折射率观测资料水平不均匀性误差的先进技术之一.利用Rueger(2002)对Smith-Weitraub方程(Smith and Weintraub，1953)的改写形式计算折射率：

N=k1P/T+k2e/T+k3e/T2,

(1)

式中N为折射率，P和e分别为大气压和水汽分压，单位为 hPa.T为大气绝对温度，单位为 K.k1、k2、k3为描述干气压项和湿气压项的常值系数.GPS 信号路径轨迹方程可用二阶微分方程表示：

(2)

其中，u=u(s)=(x(s),y(s))T是射线轨迹，n为折射指数，微分变量ds=dL/n，L是射线光学长度.方程(2)可以写成一个等价的一阶微分方程组：

du/ds=v(s),

(3a)

(3b)

其中v=du/ds，方程组的边界由u，v的初始值给出.

1.2 GNOS弯曲角正演中的水成物影响

对于中性大气，GNOS折射率进一步考虑为大气干空气分压P、温度T、水汽分压e以及液态水含量LWC(单位为g·m-3)的函数，具体表示为

N=k1P/T+k2e/T+k3e/T2+k4LWC，

(4)

k4为常值系数.根据气体状态方程P=ρdRdT及e=ρvRvT，ρd和ρv分别为干、湿空气密度，ρ为空气密度，ρd=ρ-ρv，则

N=k1Rdρ+(k′2+k3T)e/T+k4LWC

≡Ndry+Nwet+NLWC.

(5)

液态水含量LWC定义为云水和雨水含量之和(单位为g·m-3)，Rd和Rv分别为干空气和水汽气体常量，k′2=k2-k1Rd/Rv.方程(5)右侧第一项为干空气作用项(静力作用项)，第二项为湿空气作用项(非静力作用项)，第三项为液态水作用项.

饱和有云状态下，密度仅与温度及气压有关，且掩星射线在大气层中穿过的水平距离中仅有部分片断受云雨影响.定义经验参数μ表征饱和与非饱和大气的权重(也可称之为云雨区判据)，有：

N=Ndry+(1-μ)Nwet+μNsat+NLWC,

(6)

其中，Nsat=3.73×105(es/T2)，es为饱和水汽压，

es=6.112×exp[17.67(T-273.15)/(T-29.65)].

(7)

同时，饱和状态下空气密度为

ρs=P/{RdT[1+0.378es/(P-0.378es)]}-1.

(8)

通常，对于水云可采用平均统计意义的相对湿度确定参数μ，而对于冰云则可根据垂直冰水含量计算μ.但此两种方法皆不能考虑水成物的相态以及是否发生降水(雨水等)，不能包含所有水成物对折射率模拟的贡献.因此，我们直接采用所有水成物的总和作为云雨区判据，使得折射率正演计算中能够完全考虑多相态水成物的影响.因此，定义Wc为总水成物含量，即Wc=ρs(Qc+Qr+Qs+Qi+Qg)，表示单位体积内含水量总和，包括云水Qc、雨水Qr、雪水Qs、冰水Qi、霰Qg，单位皆为g·m-3，则参数μ可定义为

μ=min[1，aWc]，

(9)

此时Wc的作用仅为权重因子，已失去其作为水含量的意义，因而需进行归一化处理，a为归一化系数，同时进行强制限制，确保μ处于0到1之间.

假设大气模式云雨区中的Wc存在最大值Wcmax，则前者除以后者即可归一化. 但考虑到射线追踪算法的本质是基于非球形对称假设，在掩星剖面内自由搜索计算模拟射线路径，即弯曲角可看作射线路径上折射率的累积效应.对于对流发展深厚宽广的台风系统，当模式水平分辨率不同(所用微物理方案也随之有差异)，水成物，尤其垂直平均冰水含量的水平分布梯度(平滑度)也不同，甚至不连续，造成云雨参数对模拟折射率的贡献也不同，同化后对模式变量的调整差异较大.因此，引入模式水平分辨率参数来削弱此不确定现象.通过尝试，此处经验性做法为将模式分辨率作为Wcmax的弱约束控制参数.式(9)归一化系数定义如下：

a=1/(Wcmax+λds),

(10)

式中，λ为弱约束系数，ds为模式水平分辨率(m).

为了确定Wcmax，通过对数值预报模式中台风和锋面降水云系水成物含量总和的量级大小进行统计，发现峰值(极值)大于0.5 g·m-3的格点较少，因此取Wcmax=0.5，一方面起到归一化作用，另一方面平滑掉极端情形.同时为了确定λ，经过不同分辨率试验测试发现，取λ=0.67，基本能使得分辨率变化时增量调整大致相当.

另外，对于式(6)中的液态水作用项，有k4Lwc≥0，与式(1)对比发现，若考虑此贡献项，则云雨区中折射率模拟值恒增大，因而其取值须满足既能修正负偏差(模拟小于观测)，又不会新增正偏差(模拟大于观测)的要求.本工作据此原则进行统计计算，取值k4=1.75较为适宜，但因各模式或不同微物理过程对云雨水成物的模拟结果不尽相同，此值并不具备代表性，其他场景的实际应用中应依据背景场资料统计得出合适值.

2 FY-3c GNOS掩星弯曲角同化影响试验

2.1 FY-3c GNOS 资料

FY-3c卫星发射于2013年9月23日，是我国第二代极轨气象卫星，目标是实现全球大气和地球物理要素的全天候、多光谱和三维观测.星上搭载的GNOS提升了全球大气三维和垂直探测能力.实际检验表明，5～25 km 高度范围内,GNOS资料反演的温度误差水平优于2 K, 湿度误差优于1.5 g·kg-1(廖蜜等，2015； Liao et al., 2016).目前，GNOS掩星探测数据量仍偏少，随着北斗卫星的不断发射，以及掩星探测技术的不断发展(王先毅等，2013；乐新安等，2016)，未来GNOS掩星探测数据量有望大幅增加.在本文台风模拟试验中，为了减小资料的相关性，GNOS弯曲角资料垂直分辨率由约100 m稀疏到约200 m,并忽略300 km以内的切点飘移影响.

2.2 试验设计

为了考察观测算子改进算法的同化效果，基于WRFDA v3.7版本同化系统，利用上述算法改进后的射线追踪观测算子，发展了能够同化FY3c GNOS掩星弯曲角资料的同化模块.同时基于WRF(Weather Research and Forecasting)中尺度数值预报模式，针对2018年第24号台风“潭美”，进行GNOS弯曲角资料三维变分同化和模拟预报试验.

WRF模式模拟区域中心为(25.5°N,120.5°E),格点分布为(东西453×南北375)，水平格距15 km，垂直分层为非均匀51层，模式顶高度10 hPa.同化中，采用美国GFS(Global Forecast System)12 h预报场作为WRF区域模式背景场，同化时间窗长度为4 h，观测误差协方差利用NCEP高精度再分析资料统计模拟值与观测值之差的平方和的平均估计得到，背景误差协方差采用同化系统自身提供的基于NMC(National Meteorological Centre，美国国家气象中心)方法的统计数据，下降算法采用共轭梯度法，同化计算中数组维度约为3.4×107.另外，同化中，除了弯曲角资料，还需用到影响参数、大地水准面起伏高度、背景场温压湿及云雨水成物等参量.试验方案设计如表1所示，台风模拟初始时刻为2018年9月25日00UTC，模拟积分96 h，每3 h输出一次结果.

表1 数值试验方案Table 1 Numerical experiment schemes

分析时刻2018092500 UTC前后4 h内的GNOS弯曲角资料廓线分布如图1所示.总共10条观测廓线，每条廓线最低点资料高度差异较大，模式层内可用资料总量为1313个弯曲角观测.图中，台风标识为台风中心观测位置，圆点为资料廓线位置(颜色表征廓线起始高度)，其中2条资料廓线处于台风中心附近.选取台风中心附近右上一侧观测廓线进行分析(不妨称为“分析廓线(analyzed profile)”)，此廓线低层高度与其他廓线相比也最低.其掩星垂直剖面与FY-2f静止卫星红外1通道TBB(Temperature of Brightness Blackbody)云图的交线片段如图1b紫色线头所示，交线中点最接近台风中心，箭头表示掩星剖面从左至右的方向.可见，掩星剖面近乎穿过台风中心，且掩星无线电射线穿过台风云雨区的情况时有发生.

2.3 水成物对O-B(Observation minus Background)的影响

上述交线片段所对应掩星剖面上的模式背景场水成物WC分布如图2a所示.按照上文对水成物的定义，图中给出的水成物分布同时包含了液态水成物(云水、雨水)和固态水成物(云冰、雪、霰).图2b为“分析廓线”垂直线观测点上未考虑水成物影响的折射率模拟偏差分布(红线)和液态水含量对折射率贡献NLWC分布(绿线).折射率模拟偏差低层大于中高层，且偏差分布与NLWC分布反位相.NLWC(即K4·LWC)在11 km以下为零，向上渐增，在14 km处达到峰值，继而减小至17 km处为零.这样的分布与图2a中WC垂直分布正相关，但由于11 km以下WC主要来自于冰水含量等其他水成物而非LWC，因此相关性较弱.对于NLWC，通过系数K4的量级平衡调整作用，其与折射率模拟偏差量级相当，正负相反，故而可作为模拟值的误差缺失项予以订正而不应忽略.

继而，进一步考查水成物对射线路径上折射率模拟的影响特性.不考虑WC影响时，“分析廓线”最低点弯曲角模拟射线路径如图2a红线所示.以观测点为中心，路径两侧折射率O-B的差异(右侧减左侧)如图2c红色线所示，其在约3～10 km高度内普遍为负值.路径两侧的非对称特性(定义为射线路径同高度上右侧折射率减左侧折射率，两侧折射率存在差异则非球对称假设适用)与WC的分布相对应，说明尽管未考虑水成物影响，但由于考虑了湿气压项的贡献，一定程度上能体现大气的湿度特性，但通常GNOS信号在5 km以下的低层大气误差增大(Liao et al.，2015，2016)，导致了4 km高度层附近，非对称特性与WC的分布对应性较弱.

考虑到射线两侧水成物含量差异显著，10 km以下甚至完全晴空或完全在云雨区，因而引入射线追踪路径上的经验参数μ来增加水成物总含量的贡献，本文称之为经验云量，即式(9)所述经验参数.一旦存在水成物即表明是云雨区，则湿空气部分的贡献将分解为饱和与非饱和两部分，即(1-μ)Nwet+μNsat.增加水成物的贡献后，模拟射线路径(图2a绿线)两侧的折射率O-B差异(图2c中蓝色线)发生较大变化.一方面，缩小了4 km附近的差异，对GNOS正演折射率O-B在低层大气的改进起到一定的作用，另一方面，则增大了约6～9 km高度层内的非对称特性，这与图2b中NLWC的贡献分布更为一致.另外，尽管射线左侧存在强的WC分布，但除了约6～9 km高度上非对称性增强外，更高层差异普遍非对称性反而减弱，说明影响折射率计算结果的因素主要为以WC作为判据而实施的晴空和云雨区分离计算效应，具有复杂非线性，而非线性对应WC量值本身.

图1 初始时刻FY-3c GNOS资料 (a) 廓线水平位置分布； (b) “分析廓线”掩星剖面与FY-2f卫星红外1通道TBB云图的交线片段.Fig.1 FY-3c GNOS data of the initial time (a) Horizontal distribution of profiles; (b) Part of line that “analyzed profile” RO plane intersects to FY-2f satellite infrared 1 channel TBB cloud picture.

图2 水成物对“分析廓线”掩星观测正演计算的影响 (a) 掩星剖面片段上水成物(色阶，单位g·m-3)的分布以及两种试验方案廓线最低点资料的模拟射线路径；(b)廓线上的 模拟折射率与NLWC贡献项的对比； (c) 两种试验方案模拟射线路径两侧模拟折射率的偏差(非对称性).Fig.2 The impact of water content on forward calculation of the “analyzed profile” occultation observations (a) The distribution of water content (color filling, unit: g·m-3) on section of occultation plane and the simulated ray paths of the profile′s bottom observation in the two experiments; (b) The comparison of simulated profile refractivity with the contribution of NLWC； (c) The deviation (asymmetry) of the simulated refractivity on two sides of the simulated ray path in the two experiments.

进一步分析考虑WC影响前(图2a红线)后(图2a绿线)模拟射线路径的差异发现，因掩星剖面右侧水成物含量较小甚至晴空，考虑WC与否差异较小甚至高层几乎重合，但在云雨显著的右侧差异明显，高层受累积效应影响差异更大.总之，在考虑总含水量影响情形下，在10 km以下模拟射线路径一半处于水成物含量小于0.005 g·m-3的区域，其余部分多处云雨区域，因弯曲角是由折射率积分得到，在累积效应作用下，水成物的影响不应忽略.

图3 两个同化实验中“分析廓线”上弯曲角资料 模拟偏差随高度的变化(OMB，即O-B)Fig.3 Simulation bias of the bending angle data on the “analyzed profile” in two assimilation experiments with height (OMB, i.e. O-B)

再进一步考查折射率模拟算法改进后其效果在弯曲角正演中的累积效应体现，对图1a中“分析廓线”上所有弯曲角资料进行模拟偏差分析如图3所示.OMB表征观测与模拟的偏差(Bias),因不同高度上弯曲角数值量级差异较大，Bias取相对偏差，即(O-B)/O×100.对比未考虑水成物影响的试验方案1(VarNoWC)，引入水成物影响因子(经验参数μ)的试验方案2(VarWC)中除了8～10 km高度上的资料出现模拟偏差稍增大外，其余资料模拟偏差皆不同程度减小或保持相当，一定程度上体现了新算法的优越特性.同时发现，当高度大于10 km后，由于模拟的弯曲角射线路径上水成物为小量甚至零，偏差不再明显.某种程度上说明尽管模拟折射率的偏差未直接受WC决定，但弯曲角积分计算的累积效应却能直接体现出WC的有效性影响.

2.4 水成物对O-A(Observation minus Analysis)的影响

图4a为两个同化试验中对通过质量控制的所有GNOS弯曲角资料相对模拟偏差分层统计得出的RMSE分布.RMSE在5～7 km高度内明显减小，其余层次大致相当，且其分布与图3中的OMB分布相似.对比分析图1中参与同化的10条GNOS观测廓线，只有“分析廓线”的掩星剖面穿过台风中心附近云区，且仅有少数廓线穿越其他降水云系，因此在弯曲角的正演中考虑水成物后，所有廓线GNOS弯曲角的模拟偏差特性与“分析廓线”模拟偏差特性相似，但因平均效应，其改善程度不如图3明显.同化后模拟区域内所有参与最小化计算的弯曲角资料分析偏差OMA(相对偏差，即(O-A)/O×100)的概率密度函数(PDF)分布如图4b所示.可发现考虑水成物影响后，弯曲角OMB的PDF在[-1.5,-1.0]区间的可能性变小，与[1.0,1.5]区间的概率分布在数值上相当，成对称分布，并且最大概率的值从-0.02移动到0值，可以认为，OMA偏差较未考虑水成物影响的情形更接近于无偏的高斯分布.

2.5 GNOS资料同化对模式变量场的影响分析

图5a、5b分别是700 hPa和500 hPa上的VarNoWC试验位势高度分布，体现了台风的位置和环流分布.分析相应高度上其相对参照试验(CTRL)的增量分布(图5c、5d)发现，台风外围环境场变化较大，副热带高压整体呈减弱南退趋势，日本海东侧附近位势高度下降约13 gpm.同时，处于我国东北及中部地区的槽脊系统也存在相当量级的加强，而台风中心附近区域变化相对较小，且环境场的变化在两个高度层上体现出较好一致性，这种整层调整将更容易引导台风朝东北方向移动.进一步分析试验方案2(VarWC)试验相对于试验方案1(VarNoWC)的位势场增量(图5e、5f)，可发现，考虑了水成物影响的VarWC试验中，西太平洋副热带高压有增强西伸的调整趋势，中国大陆北部的弱高压脊后部(西北地区附近)减弱约10 gpm，而前部(东北地区附近)增强约4 gpm，即高压脊稍微向东偏移，与副高的相互作用增强.此调整趋势在700 hPa和500 hPa上皆有体现，整层的台风外围动力场调整将对台风移动路径有较大影响，相对于试验方案1，其更有利于台风减慢西伸.

参照试验初始场500 hPa上台风中心附近的温度(T)、比湿(QV)及位势高度(H)分布如图6a、6b、6c所示.可发现，背景场的分析台风已存在明显暖心结构(图6a)，台风外围至台风中心为高湿区并存在明显云系结构(图6b)，位势高度分布形成低值中心(图6c).沿图中箭头所示路径，对台风中心垂直剖面进行分析发现，分析台风本身整层暖心结构延伸到200 hPa以上(图6d)，台风中心附近高湿区域的高度也伸展至约300 hPa(图6e)，且底层至高层低压逐渐向高压转变(图6f)，体现出较好的台风天气系统结构.在此剖面上，试验方案1相对于参照试验的增量分布表明，台风中低层暖核升高(如900 hPa及500 hPa附近温度增加分别达到0.8 K)，但200 hPa高度以上的暖核结构被消弱约0.5 K(图6g).台风区域800至600 hPa高度层附近湿度增加，最大增长达到0.5 g·kg-1，且台风东侧增长范围大于西侧.台风低层东侧和中高层西侧湿度则呈现减小趋势，空间上构成明显非对称性分布(图6h).同时发现，台风中心附近400 hPa以下整层位势高度下降明显，最大下降幅度约11 gpm且位于低层， 而400 hPa往上的位势高度呈整体增高趋势，最大增加约7 gpm，可见不考虑WC影响的同化试验中台风有所增强(图6i).进一步分析试验方案2相对试验方案1的增量分布.考虑了水成物影响的试验方案2能进一步增强台风中低层暖核结构，尤其600 hPa以下整层温度升高，最大增高约0.1 K，同时200 hPa以上高度层温度普遍进一步小幅降低(图6g).湿度场的进一步调整主要出现在800至400 hPa之间(图6k)，与图6h进行比较发现，考虑水成物影响后， 800至600 hPa高度上台风中心西侧湿度的增长幅度趋缓，同时600至200 hPa之间台风中心西侧湿度的消弱程度也趋缓.位势高度场的进一步调整情况如图6l所示.与图6h相比较发现，考虑水成物影响后，台风中心低层位势高度下降使得台风稍有增加， 而200 hPa以上的位势高度增长幅度趋缓，同时台风东侧200 hPa高度以下整层位势高度的减弱程度有所缓和，这与上文图5呈现的副热带高压增强趋势是一致的.

图4 两个同化试验中所有弯曲角资料计算与观测的偏差统计 (a) 同化前相对模拟误差(O-B)的RMSE随高度分布； (b) 同化后观测与分析误差(O-A)的概率密度函数分布.Fig.4 Statistics of deviations of all calculated bending angle compared to observations in two assimilation experiments (a) Distribution of RMSE of relative simulation bias (O-B) along height before assimilation; (b) Probability density function distribution of observation and analysis error (O-A) after assimilation.

图6 初始时刻台风中心附近纬向垂直剖面(500 hPa温度、比湿、位势高度分布图(a,b,c)上箭头线所示)上的变量及增量分析：CTRL试验对应变量分布(d,e,f)；对应的VarNoWC-CTRL增量分布(g,h,i)；对应的VarWC- VarNoWC增量分布(j,k,l)(虚线表示负值，实线表示正值)Fig.6 Incremental analysis of latitudinal vertical profile along typhoon center (marked in 500 hPa temperature, specific humidity, potential height distribution (a,b,c)): CTRL experiment corresponding variable distribution (d,e,f); VarNoWC-CTRL corresponding variable incremental distribution (g,h,i); VarWC-VarNoWC corresponding variable incremental distribution (j,k,l) (dashed line indicates negative value, solid line indicates positive value)

图7为24 h预报的台风区域500 hPa垂直速度分布.参照试验(图7a)台风中心区域多为垂直静风，仅在台风中心西南位置出现明显下沉区.台风中心近外围存在较强的对流上升区，西南侧对流相对较强，达到340 cm·s-1，更外围区域垂直速度较小，接近静风.试验方案1(图7b)台风中心区域出现大范围下沉气流，强对流位置也出现调整，台风西侧对流增强，但在下沉区中存在局部上升区，中心近外围上升区对流强度也弱于参照试验，最大垂直速度仅为284 cm·s-1.进一步对试验方案2(图7c)进行分析发现，台风中心下沉区明显扩大，并消除了试验方案1中下沉区出现局部上升区的现象，较符合台风系统环流分布.另外，台风中心近外围上升区对流分布与试验方案1类似，但对流强度提升，垂直速度最大值约294 cm·s-1.试验方案2中上升对流区带再往外则出现了明显大范围下沉区，台风南北侧及西侧尤为明显.考虑到气流在台风低层辐合，环绕眼壁螺旋上升，到达流出层后分别向台风眼区和台风外围扩散并下沉，形成台风垂直环流，因此，试验方案2外围下沉区与实际情况更为吻合.

2.6 GNOS资料同化对台风模拟的影响评估

从台风模拟路径与观测实况的对比来看(图8a)，考虑了水成物影响的试验方案2，台风移动路径模拟更接近观测实况.对应的台风路径距离误差如图8b所示，在最初24 h内，试验方案2的台风路径模拟误差小于试验方案1，但两者误差皆小于25 km，改进幅度不明显.预报30～66 h内，试验方案2的台风路径误差则大于试验方案1.而72～96 h的模拟路径误差，试验方案2相对于试验方案1，有显著的改进.与参照试验相比，试验方案1和方案2的路径模拟平均距离误差分别减小了约22%和36%，即与传统GPS弯曲角射线追踪算法相比，考虑水成物影响后的新同化算法进一步取得了约14%的台风路径模拟改进效果.两个同化方案对台风中心海平面气压(图9a)的模拟结果皆稍强于参照试验，模拟过程中改善效果较稳定，且两者几乎相同.是否在GNOS弯曲角正演计算中考虑水成物影响，对海平面气压的模拟影响较小.对于台风最大风速的模拟(图9b)，两个同化试验除了51 h、72 h和90 h时次外，其余时次预报的台风最大风速皆大于参照试验，台风强度有所增强，这与海平面气压的改善效果相一致.同样地，是否考虑水成物对弯角同化的影响，对最大风速的模拟无明显的系统性影响特征.

图7 CTRL (a)、VarNoWC(b)、VarWC(c)试验中的500 hPa台风区域垂直速度24 h预报场分布 (色标两端数值分别为数据极值)Fig.7 24 h forecast of vertical velocity of typhoon area on 500 hPa in experiment CTRL (a), VarNoWC (b) and VarWC (c) (extremum values at both ends of the color scale bar, respectively)

图8 各试验方案台风模拟路径与观测实况的对比(a)以及对应的位置偏差(b)Fig.8 Comparison of typhoon simulated track and observational track of three tests (a), and the corresponding positional deviation (b)

图9 各试验方案模拟的海平面气压(a)和最大风速(b)与观测的对比Fig.9 Comparison of the simulated sea-level pressure (a) and max wind speed (b) to observation in each test

3 结论

本文在无线电掩星弯曲角正演算子中引入水成物的影响，针对台风个例，利用FY-3c GNOS弯曲角资料开展三维变分同化研究.得到结论如下：

(1)对于对流发展强盛而深厚宽广的台风天气系统，掩星路径上不同位置的垂直平均冰水含量差异显著，且存在不连续的情形，文中定义了经验参数μ表示射线路径上的经验云分布，考虑了微物理参数化和积云参数化等模式物理过程在不同模式分辨率条件下对水成物模拟结果的差异性.

(2)在GNOS弯曲角正演计算中引入掩星剖面的水成物影响因子后，影响折射率计算结果的因素主要为以水成物含量作为判据所实施的晴空和云雨区分离计算过程，而非单一对应水成物量值本身.此时考虑水成物的影响，能够对低层大气GNOS正演折射率偏差起到一定改进作用.

(3)当掩星剖面跨越一定厚度的台风区云雨大气时，水成物对掩星弯曲角正演的贡献与不引入水成物影响因子时的正演弯曲角偏差量级相当，此时水成物对掩星正演的影响不可忽略，并且O-A偏差较未考虑水成物影响的情形更接近于无偏的高斯分布.

(4)对于台风个例同化增量场的分析表明，增加水成物影响因子能进一步加强中低层台风暖中心强度，对西太平洋副热带高压、中高纬槽脊系统也出现整层趋势的调整，台风热力场和外围动力场调整对台风移动路径有较大影响.

(5)从台风模拟路径与观测实况的对比来看，考虑了水成物影响后台风移动路径更接近观测实况.与参照试验相比，忽略水成物影响和不忽略其影响，路径模拟平均距离误差分别减小了约22%和36%，即与传统GPS弯曲角射线追踪算法相比，考虑水成物影响后的新同化算法进一步取得了14%的台风路径模拟改进效果.

通过本项研究工作，初步验证了云雨天气尤其台风系统中融合水成物影响来对GNOS掩星弯曲角资料进行同化应用的有效性，为后续的理论和应用研究工作建立基础，也为数值天气预报中掩星资料的同化应用技术提供了参考.

致谢国家卫星气象中心为本研究提供了FY-3c卫星GPS掩星观测数据，本文使用的画图软件为Python语言及Matplotlib和Basemap绘图软件包，作者在此一并表示感谢．同时感谢ROM SAF(Radio Occultation Meteorology Satellite Application Facility)提供的无线电掩星数据处理包ROPP(Radio Occultation Processing Package).