王金成龔建東 鄧蓮堂

（中國氣象局?jǐn)?shù)值預(yù)報中心，北京100081）

GNSS反演資料在GRAPES＿M(jìn)eso三維變分中的應(yīng)用

王金成＊龔建東 鄧蓮堂

（中國氣象局?jǐn)?shù)值預(yù)報中心，北京100081）

為了進(jìn)一步提高GRAPES＿M(jìn)eso的分析和預(yù)報效果，該文在GRAPES＿M(jìn)eso三維變分同化系統(tǒng)中建立了同化GNSS／RO反演的大氣資料的觀測算子，實(shí)現(xiàn)了對GNSS／RO反演的大氣資料的同化應(yīng)用，并通過2013年7月1個月的同化和預(yù)報試驗(yàn)分析了GNSS／RO反演大氣資料對GRAPES＿M(jìn)eso模式系統(tǒng)分析和預(yù)報的影響。結(jié)果表明：增加了GNSS／RO反演大氣資料的同化后，GRAPES＿M(jìn)eso位勢高度場的分析誤差明顯減小，平均分析誤差減小約8%，預(yù)報誤差略有減小，平均預(yù)報誤差減小約1%；濕度場的分析誤差和預(yù)報誤差變化不明顯，常規(guī)觀測資料稀少的青藏高原地區(qū)的降水預(yù)報技巧有所提高，小雨到大雨的ETS（equitable threat score）評分提高約0.01，對全國及其他分區(qū)的降水預(yù)報技巧總體上有正效果。

GNSS掩星；反演資料；資料同化；GRAPES＿M(jìn)eso

引 言

GNSS（Global Navigation Satellite System）無線電掩星（GNSS／RO）觀測與現(xiàn)有的其他大氣觀測系統(tǒng)相比，具有觀測精度高、不受云雨等天氣影響、觀測覆蓋全球等優(yōu)點(diǎn)［1－5］，已廣泛應(yīng)用于數(shù)值天氣預(yù)報中［6－12］。GNSS／RO直接觀測及其反演的大氣資料的同化應(yīng)用對全球數(shù)值預(yù)報模式分析精度和預(yù)報技巧提 高 具 有 重 要 作 用［13－15］，近 年 我 國 學(xué) 者 在GNSS／RO資料的同化應(yīng)用方面也取得了可喜的研究成果［14－16］。目前，中國氣象局已實(shí)現(xiàn) GNSS／RO資料產(chǎn)品的業(yè)務(wù)化接收，產(chǎn)品包括：一級產(chǎn)品（L1），原始觀測的位相和振幅以及衛(wèi)星軌道和地面基站數(shù)據(jù)；二級產(chǎn)品（L2），彎角、折射率、反演的大氣氣壓和溫度及水汽等。

從理論上講，選擇何種級別的產(chǎn)品進(jìn)行同化，總原則是用最少的計(jì)算量獲得最高精度的分析結(jié)果。觀測算子簡單有利于節(jié)約計(jì)算量，同化原始觀測資料有利于提高同化精度。直接同化彎角或折射率產(chǎn)品不需要精確地衛(wèi)星軌道信息，不需要電離層模式和大氣模式層頂以上的大氣變量，其觀測算子簡單，計(jì)算量小［17］，故當(dāng)前主要業(yè)務(wù)中心均采用局地觀測算 子 同 化 彎 角 或 折 射 率 資 料［13，18－19］。 但 同 化GNSS／RO反演的氣壓、溫度和水汽資料也是一種可行的選擇。研究發(fā)現(xiàn)，GNSS／RO反演的氣壓、溫度和水汽與探空等精度較高的常規(guī)觀測資料相比具有很高的精度和一致性［13，20］，且其同化需要的觀測算子非常簡單，僅為簡單插值，計(jì)算量非常小，故對GNSS／RO反演大氣資料的同化也是值得嘗試的內(nèi)容。唐細(xì)壩等［14］實(shí)現(xiàn)了在GRAPES全球同化系統(tǒng)中同化GNSS／RO反演資料，連續(xù)循環(huán)試驗(yàn)表明：同化反演的氣壓、溫度和水汽也可以獲得較好的分析精度，并能夠顯著提高GRAPES全球模式的預(yù)報技巧。

GRAPES （Global／Regional Assimilation and PrEdiction System）是我國自主研發(fā)的數(shù)值天氣預(yù)報系統(tǒng)，已取得顯著進(jìn)步［21－31］。GRAPES區(qū)域模式（GRAPES＿M(jìn)eso）同化系統(tǒng)也已實(shí)現(xiàn)業(yè)務(wù)化，并取得了很好的業(yè)務(wù)預(yù)報效果［21］。當(dāng)前，GRAPES＿M(jìn)eso同化系統(tǒng)的研發(fā)已有較大進(jìn)展，實(shí)現(xiàn)了在模式面上的三維變分同化，預(yù)報技巧逐年提高，可同化的觀測資料包括探空、飛機(jī)報、云導(dǎo)風(fēng)、地面報、船舶報、地基GPS（Global Position System）反演的大氣可降水量和雷達(dá)風(fēng)廓線等。上述GRAPES＿M(jìn)eso三維變分同化系統(tǒng)可同化的觀測資料仍十分有限。增加GRAPES＿M(jìn)eso同化系統(tǒng)同化資料的數(shù)量是提高GRAPES＿M(jìn)eso預(yù)報技巧的重要途徑。雖然目前GRAPES全球／區(qū)域一體化的變分同化系統(tǒng)實(shí)現(xiàn)了折射率的同化，并顯著提高了GRAPES全球模式的預(yù)報技巧［15］，但GRAPES全球／區(qū)域一體化的變分同化系統(tǒng)還未應(yīng)用于GRAPES＿M(jìn)eso系統(tǒng)并實(shí)現(xiàn)業(yè)務(wù)化，特別是將其應(yīng)用于業(yè)務(wù)的GRAPES＿M(jìn)eso還有很多工作要做。因此，在GRAPES全球／區(qū)域一體化的變分同化系統(tǒng)應(yīng)用于業(yè)務(wù)GRAPES＿M(jìn)eso預(yù)報系統(tǒng)前，為了快速提高GRAPES＿M(jìn)eso同化系統(tǒng)同化資料的數(shù)量，提高區(qū)域預(yù)報系統(tǒng)分析和預(yù)報的精度，增加GNSS／RO反演的大氣資料的同化是一種快速有效方法，同時對研究GNSS無線掩星反演資料對分析和預(yù)報的影響具有重要意義。

本文在GRAPES＿M(jìn)eso同化系統(tǒng)中建立了同化GNSS／RO反演的氣壓和濕度廓線資料的觀測算子，將GNSS／RO反演的大氣氣壓和濕度廓線應(yīng)用于GRAPES＿M(jìn)eso三維變分同化系統(tǒng)中，通過數(shù)值試驗(yàn)分析同化GNSS／RO反演的大氣氣壓和濕度廓線對分析場和預(yù)報場的影響，實(shí)現(xiàn)了GNSS／RO反演資料在GRAPES＿M(jìn)eso業(yè)務(wù)同化系統(tǒng)中的應(yīng)用，提高了GRAPES＿M(jìn)eso業(yè)務(wù)預(yù)報系統(tǒng)的預(yù)報技巧。

1 模式和資料

本文采用中國氣象局自主研發(fā)并已業(yè)務(wù)化的GRAPES＿M(jìn)eso同化 預(yù)報系統(tǒng)［21－28］，模 式 為 2013年中國氣象局業(yè)務(wù)運(yùn)行版本。模式和同化系統(tǒng)的水平分辨率均為0.15°×0.15°，水平方向共502×300個格點(diǎn)，垂直方向?yàn)?2層，模式層頂約為10hPa，模式模擬區(qū)域?yàn)闁|亞區(qū)域（15°～64.35°N，70°～145.15°E）。

國家氣象信息中心業(yè)務(wù)上接收并解碼GNSS／RO觀測資料的各級產(chǎn)品，主要包括彎角、折射率和反演的大氣氣壓、溫度、濕度廓線。本文同化其中的GNSS／RO反演的大氣氣壓和濕度廓線。所有試驗(yàn)的背景場均采用美國NCEP提供的GFS全球模式6h預(yù)報場。

EC－Interim （http：／／data－portal.ecmwf.int／data／d／interim＿full＿daily）是歐洲中期天氣預(yù)報中心（ECMWF）的再分析資料，是ECMWF中期數(shù)值預(yù)報模式采用四維變分同化系統(tǒng)同化大量的觀測資料，特別是各種衛(wèi)星資料獲得的，被認(rèn)為是最接近大氣真實(shí)狀況的再分析場。因此，為了評估在同化原有資料基礎(chǔ)上加入GNSS／RO反演的大氣資料對GRAPES＿M(jìn)eso分析預(yù)報的影響，本文以EC－Interim再分析場作為參考，檢驗(yàn)GRAPES區(qū)域模式分析和預(yù)報場的精度。EC－Interim的水平分辨率為0.75°×0.75°，垂直方向?yàn)?7個等壓面層，頂層為1hPa。

2 觀測算子及相關(guān)參數(shù)設(shè)置

在已有的變分同化系統(tǒng)中，要實(shí)現(xiàn)一種新的觀測資料的同化應(yīng)用，在變分同化系統(tǒng)中需要3個方面工作：首先是建立觀測算子及其切線性算子和伴隨方程，其次是確定觀測誤差方差，另外要制定較好的質(zhì)量控制方案。

GRAPES三維變分同化系統(tǒng)采用的是增量分析方案［32］，其目標(biāo)函數(shù)表達(dá)式見文獻(xiàn)［32］。

GRAPES三維變分同化系統(tǒng)采用流函數(shù)（ψ）、勢函數(shù)（χ）、非平衡的Exner氣壓（πu）和比濕（或相對濕度）作為控制變量。

2.1 GNSS／RO反演資料的觀測算子

由于GNSS／RO反演大氣資料是以高度為垂直坐標(biāo)的大氣溫度、氣壓和濕度廓線，而GRAPES＿M(jìn)eso模式面（量綱為1的高度垂直坐標(biāo)）三維變分同化系統(tǒng)中的控制變量也是非平衡的Enxer氣壓和濕度。因此，其觀測算子僅為插值算子，包括水平插值和垂直插值。首先將模式變量從水平格點(diǎn)空間插值到觀測點(diǎn)，然后將在觀測點(diǎn)的模式變量插值到觀測的垂直層次。水平方向上采用雙線性插值方案，垂直方向采用三次樣條插值方案。

2.2 觀測誤差

新息增量法（innovation vector method，簡稱H－L方法）是 Hollingsworth等［33］提出的一種利用觀測和模式模擬之差（新息增量）統(tǒng)計(jì)背景誤差協(xié)方差和觀測誤差協(xié)方差的方法，具有同時估計(jì)背景誤差協(xié)方差和觀測誤差協(xié)方差、結(jié)果精度較高、代表性較好等優(yōu)點(diǎn)。用H－L方法估計(jì)的觀測誤差包含了測量誤差與代表性誤差，可以直接應(yīng)用于變分同化系統(tǒng)。H－L方法已被用來估計(jì)模式背景誤差協(xié)方差的三維結(jié)構(gòu)［34－35］和資料的觀測誤差［36］。

本文采用的觀測誤差方差為利用H－L方法估計(jì)COSMIC無線電掩星數(shù)據(jù)反演的大氣資料的觀測誤差方差（圖1）。

圖1 GNSS／RO反演的氣壓和相對濕度的觀測誤差Fig.1 The observation error of pressure and relative humidity retrieved by GNSS／RO data

2.3 垂直稀疏化

由于GNSS／RO反演的大氣資料的垂直分辨率為200m，因此其觀測誤差在垂直方向上的相關(guān)不能忽略；同時，模式垂直分辨率相對較低，目前GRAPES＿M(jìn)eso模式的垂直分辨率平均約為1km。由于以上兩方面原因，需要對GNSS／RO反演資料進(jìn)行垂直稀疏化，以減小觀測之間的垂直相關(guān)，并去除垂直層高頻觀測對模式的影響。研究表明，中尺度數(shù)值模式同化系統(tǒng)中在GNSS／RO資料垂直方向的稀疏化對模式預(yù)報具有一定的正貢獻(xiàn)［37］。目前，針對GNSS／RO觀測常用的稀疏化方案主要有3種：一是選取與模式層最近的觀測，二是選取模式層上下觀測的平均，三是同等間隔采樣。研究表明：第1種稀疏化方案具有較好的結(jié)果［37］，故本文采用第1種垂直稀疏化方案，即選取距離模式層最近的觀測資料進(jìn)行同化。

2.4 質(zhì)量控制方案

變分資料同化方法要求觀測資料滿足高斯分布，且其分析結(jié)果對質(zhì)量差的資料非常敏感，因此，資料的質(zhì)量控制是資料同化中的重要內(nèi)容。該文首要目標(biāo)是在業(yè)務(wù)上同化GNSS／RO反演的大氣資料，增加GRAPES＿M(jìn)eso同化觀測資料數(shù)量，提高GRAPES＿M(jìn)eso分析和預(yù)報場的精度，故暫時采用較簡單有效的質(zhì)量控制方法，即要求觀測與背景場之差（觀測增量）的絕對值小于n倍的觀測誤差方差。折射率反演的大氣廓線資料在反演和傳輸過程中引入了很多不確定性，根據(jù)質(zhì)量控制從嚴(yán)的原則，采用嚴(yán)格的質(zhì)量控制標(biāo)準(zhǔn)，設(shè)n＝2。根據(jù)經(jīng)驗(yàn)，由于GNSS／RO觀測資料在3～4km高度以下受到大氣水汽、信號的超折射等現(xiàn)象影響，精度較差，故系統(tǒng)僅同化高于地面3km以上的數(shù)據(jù)。

2.5 個例試驗(yàn)

為了驗(yàn)證觀測算子正確性并檢驗(yàn)GNSS／RO反演的大氣資料對分析結(jié)果的影響，以2013年7月12日00：00（世界時，下同）為例，只同化 GNSS／RO反演的大氣資料。2013年7月12日00：00在［－3h，3h）同化時間窗內(nèi)，模式模擬區(qū)域內(nèi)共有GNSS／RO反演的大氣廓線36條。觀測主要分布在40°N以北，特別是貝加爾湖西部地區(qū)（圖略）。

圖2是2013年7月12日00：00，同化GNSS／RO反演的大氣資料的目標(biāo)函數(shù)和其梯度范數(shù)隨迭代次數(shù)的變化曲線，表明目標(biāo)函數(shù)是收斂的，這驗(yàn)證了同化GNSS／RO反演的大氣資料觀測算子的正確性。

圖2 2013年7月12日00：00同化GNSS／RO反演的大氣氣壓和濕度的目標(biāo)函數(shù)梯度范數(shù)和目標(biāo)函數(shù)的收斂性Fig.2 Variations of the gradient norm and the costfunction for background and observation part with the iteration number for assimilating the GNSS／RO retrieved pressure and relative humidity at 0000UTC 12July 2013

2013年7月12日00：00僅同化GNSS／RO反演的大氣廓線資料個例結(jié)果表明（圖略），相比ECInterim再分析場，除了30°N以南、110°E以西區(qū)域外，背景位勢高度場偏低（偏高）的大部分區(qū)域，具有正（負(fù)）的分析增量。對于位勢高度，相比背景場，分析場更接近EC－Interim再分析場，對于比濕場也具有相同的結(jié)果。說明同化GNSS／RO反演的大氣資料的結(jié)果正確、合理，減小了模式初始場的誤差，提高了模式初始場的精度。

3 批量試驗(yàn)及結(jié)果

為了分析GNSS／RO反演的大氣資料的同化對GRAPES＿M(jìn)eso分析和預(yù)報的影響，本文進(jìn)行了1個月的對比試驗(yàn)，試驗(yàn)時間段為2013年7月1日00：00—2013年7月31日18：00，每6h進(jìn)行1次同化，均為冷啟動，背景場為NCEP GFS 6h預(yù)報場。兩組試驗(yàn)包括：①試驗(yàn)OPER：采用GRAPES＿M(jìn)eso業(yè)務(wù)系統(tǒng)，同化的資料包括探空、飛機(jī)報、地面報、船舶報、地基GPS反演的大氣可降水量和雷達(dá)風(fēng)廓線；②試驗(yàn)GNSS：在OPER基礎(chǔ)上增加了GNSS／RO反演的大氣氣壓和濕度的同化，資料截取采用6h的同化時間窗。

圖3給出了在GRAPES＿M(jìn)eso模擬區(qū)域內(nèi)2013年7月1日00：00—31日18：00業(yè)務(wù)上接收的6h同化時間窗的GNSS／RO各級產(chǎn)品廓線數(shù)量。2013年7月共接收到彎角廓線6399條，折射率廓線6174條，反演的大氣資料廓線只有3640條，平均每6h在模式區(qū)域中有約30條反演的大氣資料廓線。由圖3可以看出，反演的大氣廓線數(shù)量是彎角和折射率廓線數(shù)量的一半，明顯偏少。分析發(fā)現(xiàn)，只有COSMIC星座、GRACE（Gravity Recovery and Climate Experiment）衛(wèi)星和高分辨率雷達(dá)衛(wèi)星TerraSAR－X的無線電掩星數(shù)據(jù)在業(yè)務(wù)化系統(tǒng)提供了反演的大氣廓線，而 Metop－A／GRAS和 Metop－B／GRAS等其他衛(wèi)星的業(yè)務(wù)數(shù)據(jù)中未提供反演的大氣廓線資料。

分析試驗(yàn)GNSS中每天實(shí)際同化的GNSS／RO反演的大氣數(shù)據(jù)的數(shù)量（以標(biāo)準(zhǔn)氣壓為中心氣壓的大氣層內(nèi)總的數(shù)據(jù)量）和質(zhì)量控制剔除的資料數(shù)量（圖略），同化的質(zhì)量場（Exner氣壓）資料數(shù)量在100～300hPa最多，同時質(zhì)量控制剔除的資料量也較其他層次多。同化的濕度（相對濕度）資料數(shù)量在100～300hPa最多，經(jīng)質(zhì)量控制被剔除的相對濕度資料在100～250hPa最多，通過質(zhì)量控制剔除的資料數(shù)量達(dá)到了總資料量的1／3。2013年7月15日以后，300～700hPa資料數(shù)量明顯增多。

為了分析GNSS／RO反演的大氣資料同化在00：00，12：00（探空觀測多）和06：00，18：00（探空觀測少）分析時刻對GRAPES＿M(jìn)eso分析和預(yù)報的不同影響，分別討論了00：00，12：00和06：00，18：00分析時刻的分析和預(yù)報結(jié)果。對比圖4a和圖4b可知，GRAPES＿M(jìn)eso分析后的氣壓場比背景場更接近GNSS／RO反演的氣壓場。對比圖4c和圖4d可知，GRAPES＿M(jìn)eso分析后的相對濕度比背景場更接近GNSS／RO反演的相對濕度場。表明同化GNSS／RO反演的氣壓和濕度合理。

圖3 2013年7月1日00：00—2013年7月31日18：00業(yè)務(wù)上接收的6h同化時間窗的GRAPES＿M(jìn)eso模擬區(qū)域內(nèi)GNSS／RO資料廓線數(shù)量（a）彎角，（b）折射率，（c）反演的大氣廓線Fig.3 Number of profiles of bending angle（a），refractivity（b），retrieved atmosphere pressure，temperature and humidity（c）for［－3h，＋3h）observation window located in GRAPES＿M(jìn)eso simulated domain from 0000UTC 1July 2013to 1800UTC 31July 2013

對比圖5a和圖5b可知，GRAPES＿M(jìn)eso分析后的氣壓場比背景場更接近GNSS／RO反演的氣壓場。對比圖5c和圖5d可知，GRAPES＿M(jìn)eso分析后的相對濕度比背景場更接近GNSS／RO反演的相對濕度場。表明在06：00，18：00分析時刻，同化GNSS／RO反演的氣壓和濕度也是合理的。分別對比圖5b、圖5d和圖4b、圖4d發(fā)現(xiàn)，與00：00，12：00分析時刻相比，06：00，18：00分析時刻的分析場更接近GNSS／RO反演的大氣資料。表明GNSS／RO反演的大氣資料對06：00，18：00分析場的影響比對00：00，12：00要大，這是因?yàn)樵?6：00，18：00分析時刻幾乎沒有探空資料。

圖4 2013年7月00：00，12：00分析時刻GNSS／RO反演的氣壓及相對濕度標(biāo)準(zhǔn)差（a）反演的氣壓與背景場偏差的標(biāo)準(zhǔn)差，（b）反演的氣壓與分析場偏差的標(biāo)準(zhǔn)差，（c）反演的相對濕度與背景場偏差的標(biāo)準(zhǔn)差，（d）反演的相對濕度與分析場偏差的標(biāo)準(zhǔn)差Fig.4 The standard deviation of retrieved pressure and relative humidity from GNSS／RO at analysis time of 0000UTC and 1200UTC in July 2013（a）the standard deviation of difference between retrieved pressure and the background，（b）the standard deviation of difference between retrieved pressure and the analysis，（c）the standard deviation of difference between retrieved relative humidity and the background，（d）the standard deviation of difference between retrived relative humidity and the analysis

圖5 同圖4，但為06：00，18：00分析時刻Fig.5 The same as in Fig.4，but for analysis time of 0600UTC and 1800UTC

為了分析GNSS／RO反演的大氣資料同化對GRAPES＿M(jìn)eso模式分析場和預(yù)報場的影響，采用EC－Interim再分析資料作為參考，認(rèn)為其能代表大氣的真值，定義GRAPES＿M(jìn)eso背景場、分析場或預(yù)報場與EC－Interim再分析資料的平均距離

稱其為GRAPES＿M(jìn)eso背景場、分析場或預(yù)報場的均方根誤差，這里也簡稱為（近似）背景誤差、分析誤差和預(yù)報誤差。其中，代表在水平格點(diǎn)位置（i，j）GRAPES＿M(jìn)eso的某個大氣變量（如高度場）的背景場、分析場或預(yù)報場的值，代表 EC－Interim 再分析場相應(yīng)的值；M，N分別代表在經(jīng)向和緯向方向上GRAPES＿M(jìn)eso模擬區(qū)域內(nèi)的格點(diǎn)數(shù)目，分別為502和300。

圖6a、圖6b分別為2013年7月00：00，12：00分析時刻試驗(yàn)OPER和試驗(yàn)GNSS位勢高度場和比濕場的分析誤差。由圖6a可知，試驗(yàn)GNSS位勢高度場的分析誤差在200hPa以上均小于試驗(yàn)OPER，在200hPa以下有少部分時刻的分析誤差大于試驗(yàn)OPER。由圖6b可知，試驗(yàn)GNSS比濕場的分析誤差大約一半時刻小于試驗(yàn)OPER，另外一半時刻大于試驗(yàn)OPER的分析誤差。圖6c是位勢高度場月平均分析誤差，可以看出，試驗(yàn)GNSS的位勢高度場的月平均分析誤差明顯小于試驗(yàn)OPER；圖6d是比濕場月平均分析誤差，可以看出，試驗(yàn)GNSS比濕場的分析誤差與試驗(yàn)OPER基本相當(dāng)，沒有明顯差別。

圖6 2013年7月00：00，12：00分析時刻試驗(yàn)OPER及試驗(yàn)GNSS的分析誤差（a）位勢高度分析誤差，（b）比濕分析誤差，（c）位勢高度月平均分析誤差，（d）比濕月平均分析誤差（圖6a，6b中彩色陰影表示試驗(yàn)OPER分析誤差，等值線表示試驗(yàn)GNSS與試驗(yàn)OPER分析誤差的偏差，網(wǎng)格填充區(qū)域表示試驗(yàn)GNSS的分析誤差大于試驗(yàn)OPER）Fig.6 Variations of the analysis error for OPER and GNSS for analysis time of 0000UTC and 1200UTC in July 2013（a）analysis error of geopotential height，（b）analysis error of specific humidity，（c）monthly mean analysis error of geopotential height，（d）monthly mean analysis error of specific humidity（in Fig.6aand Fig.6b，the shaded denotes the analysis error of experiment OPER，the contour denotes the analysis error difference of GNSS to OPER，the net denotes the error of GNSS lager than that of OPER）

由圖7a可知，06：00，18：00分析時刻30～50hPa層次，試驗(yàn)GNSS位勢高度場分析誤差多大于試驗(yàn)OPER；而其他層次上，試驗(yàn)GNSS位勢高度場分析誤差僅有小部分大于試驗(yàn)OPER。由圖7b可知，100hPa以上，試驗(yàn)GNSS比濕場分析誤差均小于試驗(yàn) OPER；而100hPa以下，與00：00，12：00分析時刻的結(jié)果相似。由圖7c可知，試驗(yàn)GNSS位勢高度場月平均分析誤差除了30～50hPa以外均明顯小于試驗(yàn)OPER。對比圖7c和圖6c可知，00：00，12：00分析時刻30～300hPa層次兩組試驗(yàn)分析誤差均大于06：00，18：00，具體原因還需進(jìn)一步分析。圖7d表明，06：00，18：00分析時刻試驗(yàn)GNSS比濕場月平均分析誤差與試驗(yàn)OPER沒有明顯差別。

圖7 同圖6，但為06：00，18：00分析時刻Fig.7 The same as in Fig.6，but for analysis time of 0600UTC and 1800UTC

上述結(jié)果表明：同化GNSS／RO反演的氣壓和濕度后，明顯提高了GRAPES＿M(jìn)eso位勢高度分析場精度，而對濕度場影響不明顯?？傮w上講，增加GNSS／RO反演的大氣資料同化減小了GRAPES＿M(jìn)eso初始場誤差，有利于提高GARPES＿M(jìn)eso分析場精度。

圖8a、圖8b分別為00：00，12：00分析時刻位勢高度場24h和48h的預(yù)報誤差，圖8c、圖8d分別為位勢高度場的24h和48h預(yù)報誤差月平均值。由圖8a可知，一半以上時刻的位勢高度預(yù)報場，試驗(yàn)GNSS的預(yù)報誤差較試驗(yàn)OPER偏小。由圖8c、圖8d可知，試驗(yàn)GNSS位勢高度場24h和48h預(yù)報誤差較試驗(yàn)OPER略小。06：00，18：00分析時刻位勢高度場24h和48h預(yù)報誤差（圖略）與圖8類似，但500hPa以下，試驗(yàn)GNSS的預(yù)報誤差普遍小于試驗(yàn)OPER，相對于00：00，12：00分析時刻，06：00，18：00分析時刻試驗(yàn) GNSS在500hPa以下位勢高度場平均預(yù)報誤差減小更明顯。即同化GNSS／RO反演的大氣資料在06：00，18：00分析時刻對預(yù)報的影響更明顯。

圖8 2013年7月00：00，12：00分析時刻位勢高度場24h和48h預(yù)報誤差（a）位勢高度場24h預(yù)報誤差，（b）位勢高度場48h預(yù)報誤差，（c）位勢高度場月平均24h預(yù)報誤差，（d）位勢高度場月平均48h預(yù)報誤差（圖8a，8b中，陰影表示試驗(yàn)OPER的位勢高度的預(yù)報誤差，等值線表示試驗(yàn)GNSS與OPER預(yù)報誤差偏差，網(wǎng)格表示試驗(yàn)GNSS預(yù)報誤差大于試驗(yàn)OPER）Fig.8 Variations of 24－h(huán) and 48－h(huán) forecast errors of geopotential height for analysis time of 0000UTC and 1200UTC in July 2013（a）24－h(huán) forecast error，（b）48－h(huán) forecast error，（c）monthly mean 24－h(huán) forecast error，（d）monthly mean 48－h(huán) forecast error（in Fig.8aand Fig.8b，the shaded denotes the error of OPER，the contour denotes the forecast error difference of GNSS to OPER，the net denotes the forecast error of GNSS larger than that of OPER）

累積降水預(yù)報評分是衡量區(qū)域模式預(yù)報效果的重要指標(biāo)。圖9為00：00，12：00分析時刻全國、華南地區(qū)、長江流域、西南東部、青藏高原、華北地區(qū)、西北東部、新疆和東北地區(qū)各區(qū)域試驗(yàn)GNSS和試驗(yàn)OPER不同預(yù)報時效的24h累積降水預(yù)報的ETS（equitable threat score）評分偏差。由圖9a可知，增加GNSS／RO反演的大氣資料同化后，0～24h預(yù)報的24h累積降水ETS評分略有增加的區(qū)域和降水量級包括全國范圍的小雨和中雨，華南地區(qū)的小雨、大雨和暴雨，長江流域的中雨和大暴雨，西南東部的小雨和大雨，西北東部的小雨、中雨和大暴雨，新疆的小雨，東北地區(qū)的小雨、中雨和大雨。

圖9 2013年7月00：00，12：00分析時刻試驗(yàn)GNSS降水預(yù)報0～24h預(yù)報（a），12～36h預(yù)報（b）和24～48h預(yù)報（c）的24h累積降水ETS評分月平均值與試驗(yàn)OPER偏差Fig.9 Differences of ETS values of 0－24h（a），2－36h（b）and 24－48h（c）accumulated rainfall of GNSS to OPER for analysis time of 0000UTC and 1200UTC in July 2013 1

由圖9b可知，增加GNSS／RO反演的大氣資料同化后，12～36h的24h累積降水ETS評分增加的區(qū)域及降水量級包括全國的中雨、暴雨和大暴雨，華南區(qū)域的大雨、暴雨和大暴雨，長江流域的中雨、暴雨和大暴雨，西南東部地區(qū)的小雨和中雨，青藏高原的小雨、中雨和大雨，華北地區(qū)的中雨和大暴雨，西北東部地區(qū)的小雨、中雨、大雨、暴雨和大暴雨。由圖9c可知，增加GNSS／RO反演的大氣資料同化后，24～48h預(yù)報的24h累積降水ETS評分增加的區(qū)域和降水量級包括全國的大雨和暴雨，華南地區(qū)的小雨、大雨和暴雨，長江流域的中雨和大暴雨，西南東部地區(qū)的大雨和暴雨，青藏高原的小雨、中雨和暴雨，華北地區(qū)的大雨、暴雨和大暴雨，西北東部的大雨，東北地區(qū)的小雨、中雨、大雨和暴雨。比較圖9b和圖9c可知，增加GNSS／RO反演的大氣資料同化后，青藏高原12～36h和24～48h時刻的24h累積降水的小雨、中雨和大雨量級的ETS評分均增加，比較圖9a和圖9b可知，西北東部0～24h和12～36h時效的24h累積降水的小雨、中雨和大雨量級的ETS評分均有增加，同時12～36h預(yù)報的暴雨和大暴雨的ETS評分也有明顯增加，相對其他地區(qū)的結(jié)果增大較明顯；這可能是由于在青藏高原和西北東部試驗(yàn)OPER同化常規(guī)觀測資料量相對較少，因此在青藏和西北東部GNSS／RO反演的大氣資料同化對降水預(yù)報的影響相對明顯?？傮w上講，增加GNSS／RO反演的大氣資料同化后，青藏高原和西北東部的降水預(yù)報技巧提高明顯，新疆的降水預(yù)報技巧有所降低，而其他地區(qū)某些量級的降水的預(yù)報技巧略有提高。

圖10為2013年7月06：00，18：00分析時刻，試驗(yàn)GNSS和試驗(yàn)OPER不同預(yù)報時效的24h累積降水預(yù)報ETS評分偏差。由圖10a可知，加入GNSS／RO反演的大氣資料同化后，6～30h預(yù)報時效的24h累積降水ETS評分略有增大的區(qū)域和降水量級包括全國范圍的大暴雨，華南地區(qū)的大雨和大暴雨，長江流域的小雨和大暴雨，西南東部的小雨和暴雨，青藏高原的小雨，華北地區(qū)的大暴雨，西北東部地區(qū)的中雨，新疆的中雨，東北地區(qū)的小雨、暴雨和大暴雨。由圖10b可知，增加GNSS／RO反演的大氣資料同化后，18～42h預(yù)報時效的24h累積降水的ETS評分增大的區(qū)域及降水量級包括全國的中雨，華南區(qū)域的中雨、暴雨和大暴雨，長江流域的小雨、中雨和大雨，西南東部地區(qū)的中雨、大雨、暴雨和大暴雨，青藏高原的小雨、中雨和大雨，華北地區(qū)的暴雨和大暴雨，西北東部地區(qū)的大雨，新疆的中雨，東北地區(qū)的小雨、大雨和大暴雨。

因此，在青藏高原和西北東部等觀測資料稀少的區(qū)域，增加GNSS／RO反演的大氣資料同化后，對GRAPES＿M(jìn)eso降水預(yù)報技巧具有明顯正貢獻(xiàn)；同時對全國6～30h和12～36h預(yù)報時效的大暴雨有非常顯著的正貢獻(xiàn)。

圖10 2013年7月06：00，18：00分析時刻試驗(yàn) GNSS 6～30h預(yù)報（a），18～42h預(yù)報（b）24h累積降水ETS評分月平均值與試驗(yàn)OPER偏差Fig.10 Differences of ETS values of 6－30h（a），18－42h（b）accumulated rainfall of GNSS to OPER for analysis time of 0600UTC and 1800UTC in July 2013

4 小 結(jié)

本研究在GRAPES＿M(jìn)eso三維變分系統(tǒng)中實(shí)現(xiàn)了對GNSS／RO反演的大氣資料的同化，分析了GNSS／RO反演的大氣資料對GRAPES＿M(jìn)eso分析和預(yù)報效果的影響。結(jié)果表明：

1）業(yè)務(wù)上獲取的GNSS／RO反演的大氣資料數(shù)量明顯少于彎角和折射率，反演資料數(shù)量僅為彎角和折射率的一半。

2）增加GNSS／RO反演的大氣資料的同化顯著減小了GRAPES＿M(jìn)eso的位勢高度場分析誤差，平均分析誤差減小約8%，但對濕度場分析誤差影響不大；GRAPES＿M(jìn)eso的位勢高度場預(yù)報誤差略有減小，平均預(yù)報誤差減小約1%，濕度場的預(yù)報誤差變化不明顯。

3）2013年7月06：00，18：00分析時刻30～300hPa的位勢高度分析誤差減小量明顯大于00：00，12：00分析時刻的分析誤差，原因是06：00，18：00模式區(qū)域中幾乎沒有探空觀測。

4）GNSS／RO反演的大氣資料同化對青藏高原等常規(guī)觀測資料稀少地區(qū)的降水預(yù)報技巧的提高有明顯正貢獻(xiàn)。

本研究在GRAPES＿M(jìn)eso中實(shí)現(xiàn)了對GNSS／RO反演的大氣資料的同化應(yīng)用，改進(jìn)了GRAPES＿M(jìn)eso分析場，在一定程度上提高了GRAPES＿M(jìn)eso的預(yù)報技巧，同時為應(yīng)用GNSS／RO其他產(chǎn)品積累經(jīng)驗(yàn)，但對GNSS／RO資料同化應(yīng)用而言，折射率資料的直接同化仍然是GRAPES＿M(jìn)eso三維變分未來發(fā)展的重點(diǎn)。需要指出的是，本研究降水檢驗(yàn)所用觀測資料均為00：00—24：00的24h累積降水量，而對區(qū)域模式降水預(yù)報而言，這樣的檢驗(yàn)比較粗糙，今后工作中要進(jìn)行更精細(xì)的檢驗(yàn)。

［1］ Ware R，Rocken C，Solheim F，et al.GPS sounding of the atmosphere from low earth orbit：Preliminary results.Bull A－mer Meteor Soc，1996，77：19－40.

［2］ Wickert J，Arras C，Ao C O，et al.CHAMP，GRACE，SACC，TerraSAR－X／TanDEM－X：Science Results，Status and Future Prospects.GRAS SAF Workshop on Applications of GPSRO Measurements，2008：43－52.

［3］ Rocken C，Kuo Y H，William S S.COSMIC system description.TAO，2000，11（1）：21－52.

［4］ Anthes R A，and Coauthors.The COSMIC／FORMOSAT－3mission：Early results.Bull Amer Meteor Soc，2008，89：313－333，doi：10.1175／BAMS－89－3－313.

［5］ Kursinski E R，and Coauthors.Initial results of radio occultation observations of Earth’s atmosphere using the Global Positioning System.Science，1996，127：1107－1110.

［6］ Zou X，Kuo Y H，Guo Y R.Assimilation of atmospheric radio refractivity using a nonhydrostatic adjoint model.Mon Wea Rev，1995，123：2229－2249.

［7］ Kuo Y H，Sergey V S，Richard A A，et al.Assimilation of GPS radio occultation data for numerical weather prediction.TAO，2000，11（1）：157－186.

［8］ Liu H，Zou X，Shao H，et al.Impact of 837GPS／MET bending angle profiles on assimilation and forecasts for the period June 20－30，1995.J Geophys Res，2001，106（D23）：31771－31786.

［9］ Chen S Y，Huang C Y，Kuo Y H，et al.Assimilation of GPS refractivity from FORMOSAT－3／COSMIC using a nonlocal operator with WRF 3DVAR and its impact on the prediction of a typhoon event.Terr Atmos Ocean Sci，2009，20（1）：133－154.

［10］ Ma Z Z，Kuo Y H，Wang B，et al.Comparison of local and nonlocal observation operators for the assimilation of GPS RO data with the NCEP GSI System：An OSSE study.Mon Wea Rev，2009，137：3575－3587，doi：0.1175／2009MWR2809.1.

［11］ Huang C Y，Kuo Y H，Chen S Y，et al.Impact of GPS radio occultation data assimilation on regional weather predictions.GPS Solut，2010，14：35－49，doi：10.1007／s10291－009－0144－1.

［12］ Ma Z Z，Kuo Y H，Martin F R，et al.Assimilation of GPS radio occultation data for an intense atmospheric river with the NCEP regional GSI system.Mon Wea Rev，2011，139：2170－2183，doi：10.1175／2011MWR3342.1.

［13］ Cucurull L，Derber J C.Operational implementation of COSMIC observations into NCEP’s Global Data Assimilation System.Wea Forecasting，2008，23：702－711.

［14］ 唐細(xì)壩，薛紀(jì)善.COSMIC資料在GRAPES全球三維變分同化系統(tǒng)的初步研究.熱帶氣象學(xué)報，2009，25（5）：521－531.

［15］ Liu Y，Xue J S.Operational implementation of the assimilation of global navigation satellite radio occultation observations into GRAPES data assimilation system.J Meteor Res，2014，28，doi：10.1007／s13351－014－4028－0.

［16］ 趙德顯，郁紅弟，沈桐立.GPS折射率資料在梅雨鋒暴雨數(shù)值模擬中的應(yīng)用.氣象，2011，39（12）：1511－1518.

［17］ Syndergaard S，Kuo Y H，Lohmann M S.Atmosphere and Climate，Chapter：Observation Operators for the Assimilation of Occultation Data into Atmospheric Models：A Review.Berlin：Springer Berlin Heidelberg，2006：205－224.

［18］ Poli P，Healy S B，Dee D P.Assimilation of Global Positioning System radio occultation data in the ECMWF ERA－Interim reanalysis.Q J R Meteorol Soc，2010，136：1972－1990.

［19］ Rennie M P.The impact of GPS Radio occulation assimilationat the Met Office.Q J R Meteorol Soc，2010，136：116－131.

［20］ Wang B R，Liu X Y，Wang J K.Asessment of COSMIC radio occultation retrieval product using global radiosonde data.Atmos Meas Tech，2013，6：1073－1083，doi：10.5194／amt－6－1073－2013.

［21］ 薛紀(jì)善，陳德輝.數(shù)值預(yù)報系統(tǒng)GRAPES的科學(xué)設(shè)計(jì)與應(yīng)用.北京：科學(xué)出版社，2008.

［22］ 陳德輝，楊學(xué)勝，胡江林，等.多尺度通用動力模式框架的設(shè)計(jì)策略.應(yīng)用氣象學(xué)報，2003，14（4）：452－461.

［23］ 黃麗萍，伍湘君，金之雁.GRAPES模式標(biāo)準(zhǔn)初始化方案設(shè)計(jì)與實(shí)現(xiàn).應(yīng)用氣象學(xué)報，2005，16（3）：374－384.

［24］ 萬齊林，薛紀(jì)善，陳子通，等.雷達(dá)TREC風(fēng)的三維變分同化應(yīng)用與試驗(yàn).熱帶氣象學(xué)報，2005，21（5）：449－457.

［25］ 伍湘君，金之雁，黃麗萍，等.GRAPES模式軟件框架與實(shí)現(xiàn).應(yīng)用氣象學(xué)報，2005，16（4）：539－546.

［26］ 胡江林，沈?qū)W順，張紅亮，等.GRAPES模式動力框架的長期積分特征.應(yīng)用氣象學(xué)報，2007，18（3）：276－284.

［27］ 伍湘君，金之雁，陳德輝，等.新一代數(shù)值預(yù)報模式GRAPES的并行計(jì)算方案設(shè)計(jì)與實(shí)現(xiàn).計(jì)算機(jī)研究與發(fā)展，2007，44（3）：510－515.

［28］ 萬齊林，薛紀(jì)善.曲率修正線性平衡方程及其在變分同化風(fēng)壓約束中的應(yīng)用.熱帶氣象學(xué)報，2007，23（5）：417－423.

［29］ 王莉莉，陳德輝，趙琳娜.GRAPES氣象－水文模式在一次洪水預(yù)報中的應(yīng)用.應(yīng)用氣象學(xué)報，2012，23（3）：274－284.

［30］ 王瑞春，龔建東，張林.GRAPES變分同化系統(tǒng)中動力平衡約束的統(tǒng)計(jì)求解.應(yīng)用氣象學(xué)報，2012，23（2）：129－138.

［31］ 王雨，李莉.GRAPES＿M(jìn)eso V3.0模式預(yù)報效果檢驗(yàn).應(yīng)用氣象學(xué)報，2010，21（5）：524－534.

［32］ 莊世宇，薛紀(jì)善，朱國富，等.GRAPES全球三維變分同化系統(tǒng)－基本設(shè)計(jì)方案與理想試驗(yàn).大氣科學(xué)，2005，29（6）：872－884.

［33］ Hollingsworth A，Lnnberg P.The statistical structure of shortrange forecast errors as determined from radiosonde data.Part I：The wind field.Tellus，1986，38：111－136.

［34］ Xu Q，Li W，Andrew V T，et al.Estimation of three－dimensional error covariances.Part I：Analysis of height innovation vectors.Mon Wea Rev，2001，129：2126－2135.

［35］ 龔建東，魏麗，陶士偉，等.全球資料同化中誤差協(xié)方差三維結(jié)構(gòu)的準(zhǔn)確估計(jì)與應(yīng)用I：觀測空間協(xié)方差的準(zhǔn)確估計(jì).氣象學(xué)報，2006，64（6）：669－684.

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

［37］ Hsieh M E，Chang L Y，Hsiao L F，et al.Impact of Quality Control and Data Thinning of GPS RO Data in WRF－Var on Typhoon Track Forecast.Sixth FORMOSAT－3／COSMIC Data Users’Workshop，NCAR Center Green Campus，2012.

Operational Assimilation of Data Retrieved by GNSS Observations into GRAPES＿M(jìn)eso 3DVar System

Wang Jincheng Gong Jiandong Deng Liantang
（Numerical Weather Prediction Center of CMA，Beijing100081）

Radio occultation（RO）observations using the Global Navigation Satellite System （GNSS）provides valuable data to support operational numerical weather prediction，and it is proved that the assimilation of RO data has the potential to significantly improve the accuracy of global and regional meteorological analysis and weather prediction.RO observations have many advantages such as high precision，global coverage and high vertical resolution compared with other observations.RO observations provide phase（Level 1a），bending angle（Level 1b），refractivity（Level 2a），retrieved pressure，temperature and humidity profiles（Level 2b），all of which can be assimilated into numerical model.Bending angle and refractivity are proved better choices for assimilation as observation operators are less complicated，and have no disadvantages associated with the modeling of ionospheric effects in the assimilation model.However，assimilation of water vapor，temperature or pressure derived from RO observations have advantages that data forms are same as model variables，and it is proved that assimilation of retrieved data also improve the accuracy of global analyses and forecasts significantly.

In the operational GRAPES （Global／Regional Assimilation and Prediction System）regional（GRAPES＿M(jìn)eso）3－dimensional variational（3DVar）system，the analysis is performed in model vertical level，but it can only assimilate the upper－level wind，pressure and humidity from radiosonde report（TEMP），upper－level wind and humidity by AIREP，cloud drift wind（SATOB），the pressure from surface station report over land（SYNOP）and over sea，integrated column precipitable water vapor（IPW）retrieved by ground－based GPS observations and radar wind profiles.In order to assimilate RO retrieved atmosphere data and address their impacts on analyses and forecasts，observation operators of retrieved humidity and pressur eassimilation are developed first，and then quality control and vertical thinning scheme are formed and discussed.Results of one month experiments show that the assimilation of GNSS／RO retrieved pressure and humidity data have considerable positive impacts on analyses of geopotential height，minor positive impacts on humidity analyses，geopotential height，humidity forecasts，rainfall forecasts，and major positive impacts on rainfall forecast for Qinghai－Tibet region.

GNSS radio occultation；retrieved data；data assimilation；GRAPES＿M(jìn)eso

王金成，龔建東，鄧蓮堂.GNSS反演資料在GRAPES＿M(jìn)eso三維變分中的應(yīng)用.應(yīng)用氣象學(xué)報，2014，25（6）：654－668.

2014－05－13收到，2014－09－17收到再改稿。

國家國際科技合作專項(xiàng)（2011DFG23210），公益性行業(yè)（氣象）科研專項(xiàng)（GYHY201406011），國家氣象中心青年基金項(xiàng)目“GNSS無線電掩星反演資料在GRAPES區(qū)域三維變分中的同化應(yīng)用”

＊email：wangjc＠cma.gov.cn