GPS、洋面散射计等资料质量控制

观测误差统计

在变分同化中，背景误差和观测误差共同决定了背景场和观测资料在分析场中所占的比例。观测误差的合理给定是资料同化的重要环节，对GPS ZTD资料同化同样如此。Vedel H等（2004）在HIRLAM模式中采用1cm作为地中海西部GPS ZTD资料的观测误差并进行同化。仲跻芹（2017）则采用1.25cm作为华北地区GPS ZTD资料的观测误差进行同化应用。然而不同地区资料质量存在差异，质量控制对资料观测误差也有一定的影响。本章采用多种观测误差统计方法对于质量控制后的GPS ZTD资料进行观测误差统计，并比较分析不同观测误差统计结果。

Hollingsworth方法

方法介绍

Hollingsworth方法一开始被应用于估计无线电探空仪网络的预报误差（Hollingsworth等，1986）。其主要思路是：首先，基于真实背景误差在空间上相关，而观测误差在空间上是不相关的假设。通过计算预报场和观测离差协方差的相关程度作为距离的函数来估计观测误差。通过将预报场和观测离差协方差的相关性从非零分离外推到零分离处的背景误差和观测误差来估计观察误差，使得零分离处的预报场和观测离差方差被分成空间相关分量和空间不相关分量，假设后者为观测误差部分。该方法同时还需假设观测误差和背景误差是不相关的。

该方法的变体也用于通过大气运动矢量AMV和无线电探空仪之间的差异来估计大气运动矢量（AMV）的空间误差相关性（Bormann等，2003）。Bormann N 和Bauer P（2010）将该方法作为误差统计方法之一估计了AMSU-A微波等卫星资料的观测误差。该方法假设观测误差空间不相关，但由于卫星资料空间分布比较密集，往往需要考虑观测误差的空间相关性，因此他们指出该方法在估计卫星观测误差时存在一定局限性。

此外当研究区域的数据比较密集时，与背景误差相关性的长度尺度相比最短间隔距离往往相当小。因此，可以不采用这种相关性函数进行外推，而是在第一个具有足够观测矢量最小距离的非零点中分离观测误差和背景误差。由于预报离差协方差倾向于随着距离的减少而增加，忽略相关性函数的使用可能会导致预报场离差方差的空间相关部分的低估，因此高估了观测误差。

3.1.2观测误差统计分析

在采用Hollingsworth方法统计GPS ZTD资料的观测误差时，首先计算GPS ZTD资料相对于3小时预报场的离差，然后每隔100km计算不同距离上离差的相关性。最后通过引入拟合的函数将相关性从非零外推到零距离，基于真实背景误差在空间上相关而观察误差在空间上是不相关的假设，对观测误差和背景误差进行分离。图3.1显示了不同距离上离差的相关性及不同间隔距离上匹配的离差数量。从图3.1（a）可以看到在400km内基本上随着距离的增加相关性呈现减小的趋势，400km以外相关性变化不大。从图3.1（b）可以看在不同间隔距离上匹配到的离差数量大部分都在1000个以上，这基本保证了在在统计不同间隔距离相关性的稳定性。

本研究通过引入拟合函数外推得到在零距离处的相关性为0.544，而且计算得到GPS ZTD观测和3小时预报场的离差的均方根为1.698cm，因此最终统计得到的GPS ZTD观测误差为0.78cm。

（a）不同距离上离差相关性及（b）不同间隔距离上匹配的离差数量

Desroziers diagnostic（2001）

方法介绍

Desroziers和Ivanov（2001）提出了一种基于最小成本函数的最优性标准来调整观测误差协方差的比例系数的方法，该方法假设相关性在最初假定的观测误差协方差中准确表示。Chapnik 等（2004）阐述了该方法的一些性质，并展示了该方法用于卫星观测的初步结果，处理ARPEGE同化系统中实际卫星辐射的试验表明，调谐系数与标准偏差的实际值相关。Chapnick等（2006）提出该方法的新实现方式，其中包括用于大矩阵迹线的新随机化方案，并研究该方法调整的方差在同化系统中的一些初步结果。

对于代价函数 \(J_j(x^a)\) 而言其期望值：

\[\begin{split}E(J_j) & = \frac{1}{2} E\{(\Gamma_j x^a - y_j^o)^\intercal R_j^{-1} (\Gamma_j x^a - y_j^o)\} \\ & = \frac{1}{2} E\{Tr ( R_j^{-1} (\Gamma_j x^a - y_j^o)^\intercal (\Gamma_j x^a - y_j^o))\} \\ & = \frac{1}{2} E\{Tr ( R_j^{-1} (\Gamma_j \varepsilon^a - \varepsilon_j^o)^\intercal (\Gamma_j \varepsilon^a - \varepsilon_j^o))\} \\ & = \frac{1}{2} Tr ( R_j^{-1} [\Gamma_j E(\varepsilon^a \varepsilon_j^{a\intercal}) \Gamma_j^\intercal + E(\varepsilon^a \varepsilon_j^o)^\intercal - E\{(\Gamma_j \varepsilon^a) \varepsilon_j^{o\intercal}\} - E\{(\Gamma_j \varepsilon^a)^\intercal \varepsilon_j^o)\}])\end{split}\]

其中为 \(x^a\) 分析场，为 \(y_i^o\) 观测，\(R_j\) 为观测算子，\(\Gamma_j\) 为线性观测算子，\(\Gamma_j x^a - y_j^o = \Gamma_j x^a + \Gamma_j x^t -\Gamma_j x^t - y_j^o = \Gamma_j \varepsilon^a - \varepsilon_j^o\) ,且 \(\{E(Tr(.)\} = Tr\{E(.)\}\) 。由于 \(\varepsilon^a = \varepsilon^b + \mathrm{K} \{H(x^t) + \varepsilon^o - H(x^b)\} = \mathrm{K} \varepsilon^o\)，其中 \(\mathrm{K}\) 为增益矩阵，因此可得 \(E(\varepsilon^a \varepsilon^{o\intercal}) = \mathrm{K} E(\varepsilon^o \varepsilon^{o\intercal}) = \mathrm{P}^a H^\intercal \; , \; E\{\Gamma_j \varepsilon^a) \varepsilon_j^{o\intercal}\} = \Gamma_j E(\varepsilon^a \varepsilon^{o\intercal}) \mathrm{P}_j^\intercal = \Gamma_j \mathrm{P}^a H^\intercal \mathrm{P}_j^\intercal = \Gamma_j \mathrm{P}^a \Gamma_j^\intercal\) , \(\mathrm{P}_j\) 为投射算子，\(\varepsilon_j^o = \mathrm{P}_j \varepsilon^o\) ，\(\mathrm{P}^a\) 代表由整组观察结果分析得到的估计分析的误差协方差矩阵。\(E(\varepsilon^a \varepsilon_j^o)^\intercal = R_j\) ，所以有：

\[E(J_j) = \frac{1}{2} Tr \{ R_j^{-1} (R_j - \Gamma_j \mathrm{P}^a \Gamma_j^\intercal)\} = \frac{1}{2} \{m_j - Tr(R_j^{-1} \Gamma_j \mathrm{P}^a \Gamma_j^\intercal)\}\]

其中 \(\mathrm{P}^a\) 表示分析误差协方差，\(\Gamma_j\) 表示线性观测算子，\(R_j\) 表示为观测误差协方差，\(Tr\) 表示为矩阵的迹。所以：

\[E(J^o) = \frac{1}{2} \{ p - Tr(R^{-1} H \mathrm{P}^a H^\intercal) \} = \frac{1}{2} \{ p - Tr(\mathrm{K} H) \}\]

其中 \(\mathrm{K} = R^{-1} H^\intercal \mathrm{P}^a\) , \(J^o\) 为观测部分的代价函数。

\[J(\delta x) = \frac{1}{(s^b)^2} J^b (\delta x) + \frac{1}{(s^o)^2} J^o (\delta x)\]

其中, 可以得到：

\[\frac{1}{(s^o)^2} J^o (\delta x) = E(J^o) = \frac{1}{2} \{ p - Tr(\mathrm{K} H) \}\]

最终可以得到

\[(s^o)^2 = \frac{2 J^o}{\{p - Tr(\mathrm{K} H)\}}\]

\(p\) 代表了观测集的观测数，从这里可以算出合适的观测误差，\(Tr(HK)\) 则需要通过randomized方法进行估计：

(1)\[Rand(Tr(KH)) = Rand( Tr(R^{-1/2} K H R^{1/2}) ) = (R^{-1/2} \zeta )^\intercal K H R^{1/2} R^{-1/2} \zeta\]

(2)\[H \delta x_{y^o + \delta y^o}^a - H \delta x_{y^o}^a \; ; \; H K \delta y^o\]

(3)\[\delta y^o = R^{-1/2} \zeta\]

其中 \(\zeta\) 是具有标准高斯分布（均值为0和方差为1）的 \(p\) 维向量，\(\delta x_{y^o + \delta y^o}^a\) 和 \(\delta x_{y^o}^a\) 分别为对观测进行扰动后的分析增量和不对观测进行扰动后的分析增量。

由 (1) ，(2) 和 (3) 可以得到：

\[Rand(Tr(KH)) = (R^{-1/2} \zeta)^\intercal (H \delta x_{y^o + R^{1/2} \zeta}^a - H \delta x_{y^o}^a)\]

观测误差统计分析

本研究通过Desroziers（2001）来统计观测误差时，发现对于不同的极小化化过程，由于同化中观测数量及质量的差异，每次观测误差统计结果存在差异。图3.2显示了不同同化过程的观测误差情况，可以看到绝大多数观测误差结果位于在0.5-0.8之间，数值波动较小；而较少数量的结果波动较大，数值大于1cm。基于上述情况，个别时次的结果甚至超过了3cm，本研究将所有的结果的均值作为最终的观测误差结果，最终计算得到GPS ZTD观测误差为0.84cm。

不同同化过程的观测误差结果

Desroziers diagnostic（2005）

方法介绍

Desroziers等(2005)提出了通过OMB和OMA诊断观测误差的方法，并表明即使初始假设的观测误差是不相关的，只要真正的背景误差和真实的观测误差具有足够不同的相关结构，该方法具有检索观测误差的空间相关结构的能力。该方法已被用于估计观察误差和通道间误差相关性（M’enard等，2009; Stewart等，2009）。Bormann N 和Bauer P（2010）将该方法作为误差统计方法之一估计了AMSU-A微波等卫星资料的观测误差。Wei Han等（2007），在Desroziers (2005)的理论基础上实现了新的误差算法调整，该方法在变分同化的背景和观察误差的实际调整中被证明是可行的，并被应用于中国气象局（CMA）新一代NWP系统GRAPES中。Cordoba M等（2017）年将该方法用于云导风资料观测误差估计，并讨论了云导风观测误差的相关性。

理论及公式推导

\[x^a = x^b + \delta x^a = x^b + K d_b^o\]

其中 \(x^b\) 为背景场，\(\delta x^a\) 为分析增量，\(d_b^o\) 为观测和背景场对应部分差异，\(K\) 为增益矩阵。

\[K = \mathrm{B} H^\intercal ( H \mathrm{B} H^\intercal + \mathrm{R} )^{-1}\]

\(H\) 为非线性观测算子，\(\mathrm{R}\) 为观测误差，\(\mathrm{B}\) 为背景误差协方差。

(4)\[d_b^o = y^o - H(x^b) = y^o - H(x^t) + H(x^b) -H(x^t) \; ; \; e^o - \mathrm{H} e^b\]

\[d_a^o = y^o - H(x^b + \delta x^a) \; ; \; y^o - H(x^b) - \mathrm{H}Kd_b^o = \mathrm{R}(\mathrm{H} \mathrm{B} \mathrm{H}^\intercal + \mathrm{R})^{-1} d_b^o\]

其中 \(x^t\) 为未知的大气真实状态，\(e^o\) 为观测误差向量，\(e^b\) 为背景误差向量。

采用线性统计期望算子:math:E 并假设观测向量和背景误差向量不相关时，由:eq:f313 可以得到：

(5)\[E[d_b^o(d_b^o)^\intercal] = E[e^o(e^o)^\intercal] + \mathrm{H} E[e^b(e^b)^\intercal] \mathrm{H}^\intercal\]

(6)\[E[d_b^o(d_b^o)^\intercal] = \mathrm{H}\mathrm{B}\mathrm{H}^\intercal + \mathrm{R}\]

采用线性统计期望算子 \(E\) 并假设观测向量和背景误差向量不相关时，由:eq:f315 和 (6) 可以得到

\[E[d_b^o(d_b^o)^\intercal] = \mathrm{R} ( \mathrm{H}\mathrm{B}\mathrm{H}^\intercal + \mathrm{R} )^{-1} E[d_b^o(d_b^o)^\intercal]\]

\[E[d_b^o(d_b^o)^\intercal] = \mathrm{R}\]

所以观测误差可以通过OMB和OMA来进行统计分析。

观测误差统计分析

本研究采用WRF的三小时预报作为背景场，并计算背景场与GPS ZTD观测之间的OMB和OMA来统计指控后GPS ZTD资料观测误差。本研究统计2018年7月份一个月的GPS ZTD的资料，如表格4.1所示，在7月份参与统计的GPS ZTD观测资料的数量一共有21974个，通过Desroziers (2005)方法统计的观测误差为0.779cm。

表3.1 Desroziers (2005)方法观测误差情况

观测数量	观测误差（cm）
21974	0.779

月站点最大OMB标准差方法

该采用当月站点中OMB标准差最大值作为GPS ZTD资料观测误差，是一种相对比较保守的观测误差方法，该方法可以通过质控过程中的结果中得到，在质控过程中可以得到afqc4meanbzc.dat文件。图4.3 进过质控后不同站点7月份OMB的标准差情况，剔除超过红线的站点后可以得到对的标准差为1.99cm。

不同站点7月份OMB的标准差情况

本章小结

本章节分别通过WRF-3DVar中的Hollisworth （1986）方法， Desroziers G （2001）方法、Desroziers G （2005）方法以及月站点最大OMB标准差方法对GSPS ZTD的观测误差进行统计。

表3.2显示了用不同的观测误差统计方法统计GPSZTD资料观测误差结果，可以看到不同方法的结果比较相近，Hollisworth方法统计得到的观测误差为0.78cm，Desroziers G （2001）方法统计得到的观测误差为0.84cm，Desroziers G （2005）方法统计得到的观测误差为0.779cm，而月站点最大OMB标准差方法得到的观测误差1.99cm。

表3.2 不同方法统计GPS ZTD观测误差结果

方法名称	观测误差（cm）
Hollisworth	0.78
Desroziers G （2001）	0.84
Desroziers G （2005）	0.779
最大标准差	1.99

同化及预报试验分析

为了检验统计的GPS ZTD资料观测误差统计结果对该资料的同化效果以及对短期预报的影响，本章开展了2018年夏季5次降水个例的循环同化及预报试验。从前一章可以看出，前三种方法得到的观测误差结果比较相近，所以本研究选定前三种方法的均值0.8cm和月站点最大OMB标准差方法的1.99cm，分别作为GPS ZTD资料的观测误差进行同化测试，并比较分析两种观测误差以及是否同化GPS ZTD对降水预报的影响。

同化及预报试验设计

所有试验都使用WRF 3.9.1版本（Skamarock et al., 2008）。试验中使用了两层双向嵌套网格。外网格的分辨率为9km，共701×601个模式格点；内网格的分辨率为3km，共691×703个格点（图4.1）。垂直方向内外网格都使用了50层，模式层顶为10 hPa。同时物理参数化方案选取了：WSM6微物理参数化方案(Hong and Lim, 2006)；RRTM长波辐射方案(Mlaweret al. 1997); Dudhia 短波辐射方案(Dudhia, 1989)； YSU行星 边界层方案 (Hong et al., 2006); Noah陆面方案(Tewari et al., 2004).

图4.1 试验区域图

在循环同化试验中，每三个小时产生一个分析场。3h和24h预报所需的边界条件都来自NCEP的全球预报系统的分析场。下一时刻的背景场则为上一时刻的3h预报场。静止背景误差协方差则利用前一个月WRF的12和24h预报差通过NMC算法统计得到。

本章一共设计了几组同化试验来评估同化GPS ZTD资料对本次强降水过程的预报和分析的影响。第一组试验为控制试验（简称CTRL），主要同化了来自卫星反演云导风、探空、船舶、飞机报以及海洋和陆地上的各测站资料（图4.2）。第二组试验为GPS ZTD同化试验再GPS ZTD资料观测误差采用0.8cm作为观测误差（简称GPS08），该试验所同化的常规资料与CTRL试验相同，只是在外网格上增加了GPS ZTD资料同化，第三组试验为GPS ZTD同化试验再GPS ZTD资料观测误差采用1.99cm作为观测误差（简称GPS1.99），其他试验设置和GPS08试验相同。

本章开展了2018年夏季5次降水个例的循环同化及预报的对比试验，每个个例分别从预报起点前18小时开始进过6小时的spin up，然后经过5次3小时连续循环同化，然后进行24小时预报。

图4.2 试验所同化得到常规观测资料

降水预报比较分析

（1）个例1

图4.3为2018年8月7日00时到8日00时的24小时累计降水情况（a为观测降水，b为CTRL试验预报降水，c为GPS08试验预报降水，d为GPS1.99试验预报降水）。可以看到在这次个例中在采用不同的观测误差同化GPS ZTD资料对于降水会产生一定的影响，减少了辽宁附近的降水的空报和过包，其中采用0.8cm观测误差试验降低辽宁附近的过报更加明显。

为2018年8月7日00时到8日00时的24小时累计降水

图4.4为2018年8月7日00时的可降水量，a为CTRL试验，b为GPS08试验-CTRL试验,c为GPS1.99试验-CTRL试验。从图中看到同化GPS ZTD资料后的差异主要集中辽宁附近，这是应为辽宁附近的资料是有辽宁省局单独解码得到资料质量相对好，通过质控的观测数量较多。可以看到GPS08试验和GPS1.99试验在40°N，118°E附近的可降水明显比CTRL试验低，这样导致GPS08试验和GPS1.99试验的24小时降水预报在这附近的空报和过报

为2018年8月7日00时的可降水量，a为CTRL试验，b为GPS08试验-CTRL试验,c为GPS1.99试验-CTRL试验

（2）个例2

图4.5为2018年8月14日00时到15日00时的24小时累计降水情况（a为观测降水，b为CTRL试验预报降水，c为GPS08试验预报降水，d为GPS1.99试验预报降水）。可以看到在这次个例中在采用不同的观测误差同化GPS ZTD资料对于降水会产生一定的影响，减少了辽天津附近降水的过报。同时，三组试验对渤海地区的降水强度都偏弱，位置偏北，但是GPS08试验对渤海的降水预报在降水强度和降水位置上和实况更加接近。

图4.5 为2018年8月14日00时到15日00时的24小时累计降水

图4.6为2018年8月14日00时的可降水量，a为CTRL试验，b为GPS08试验-CTRL试验,c为GPS1.99试验-CTRL试验。可以看到在同化GPSZTD资料后对于38°N，116°E附近的可降水量存在较明显的负增量，得益于可降水的减少GPS08试验和GPS1.99试验在该地附近的降水与的过报得到一定的改善。

图4.6为2018年8月14日00时的可降水量，a为CTRL试验，b为GPS08试验-CTRL试验,c为GPS1.99试验-CTRL试验

（3）个例3

图4.7为2018年8月20日00时到21日00时的24小时累计降水情况（a为观测降水，b为CTRL试验预报降水，c为GPS08试验预报降水，d为GPS1.99试验预报降水）。可以看到在这次个例中在采用不同的观测误差同化各组试验基本都能模拟出降水的基本形式，三组试验都存在较明显的过报。GPS08试验对150mm以上的强降水的范围更小，和实况更加接近，减弱的过报。

图4.7为2018年8月20日00时到21日00时的24小时累计降水

图4.8为2018年8月20日00时的可降水量，a为CTRL试验，b为GPS08试验-CTRL试验,c为GPS1.99试验-CTRL试验。可以看到同化GPSZTD后对可降水量影响主要集中在辽宁附近以及北京附近，主要是这些地区资料数量较多。在本次个例中降水量和预报的24小时降水并没有很好对应，可能是由于24小时预报时间过长导致。

图4.8 为2018年8月20日00时的可降水量，a为CTRL试验，b为GPS08试验-CTRL试验,c为GPS1.99试验-CTRL试验

（4）个例4

图4.9为2018年7月12日12时到13日12时的24小时累计降水情况（a为观测降水，b为CTRL试验预报降水，c为GPS08试验预报降水，d为GPS1.99试验预报降水）。可以看到在这次个例中在采用不同的观测误差同化各组试验基本都能模拟出降水的基本形式，三组试验都存在较明显的过报。GPS08试验在渤海地区降水预报的范围和强度比CTRL试验和GPS1.99试验都更好。同时对于北京附近的降水，CTRL试验预报降水位置偏西，GPS08试验和GPS1.99试验预报位置更好，但是在预报强度偏弱。

图4.9 为2018年7月12日12时到13日12时的24小时累计降水

图4.10为2018年7月12日12时的可降水量，a为CTRL试验，b为GPS08试验-CTRL试验,c为GPS1.99试验-CTRL试验。可以看到该次个例中可降水量的高值和降水区域对应较好。在同化GPS ZTD后GPS08试验相对于CTRL试验在40°N，30°E附近可降水量由明显的负增量，而GPS1.99试验在该地区也有负增量但是没有GPS08试验明显。所以在该地区GPS08试验相对于CTRL试验减弱了空报，GPS1.99试验则介于GPS08试验和CTRL试验之间。

图4.9 为2018年7月12日12时的可降水量，a为CTRL试验，b为GPS08试验-CTRL试验,c为GPS1.99试验-CTRL试验

（5）个例5

图4.10为2018年7月10日21时到11日21时的24小时累计降水情况（a为观测降水，b为CTRL试验预报降水，c为GPS08试验预报降水，d为GPS1.99试验预报降水）。GPS08试验在渤海地区降水预报的范围和强度比CTRL试验和GPS1.99试验都更好。CTRL试验和GPS1.99试验在长春地区存在明显空报，而GPS08试验对于长春的降水过报有明显改善。三组试验对于辽宁附近都存在过报但是GPS08降水范围和实况最为接近。

图4.10 为2018年7月10日21时到11日21时的24小时累计降水

图4.11为2018年7月10日21时的可降水量，a为CTRL试验，b为GPS08试验-CTRL试验,c为GPS1.99试验-CTRL试验。GPS08试验相对于CTRL试验的可降水的差异明显改与GPS1.99试验，主要是应为GPS08试验同化GPSZTD资料时观测误差较小。可以看到在44°N，126°E附近GPS0试验的可降水量低于CTRL试验，所以在该地区GPS08试验的空报明显小于CTRL试验。

图4.11 为2018年7月10日21时的可降水量，a为CTRL试验，b为GPS08试验-CTRL试验,c为GPS1.99试验-CTRL试验

图4.12为5个个例的降水的平均FSS评分，可以看到在所有阈值GPS08试验的FSS评分都是最高，GPS1.9试验和CTRL试验比较接近且略低于CTRL试验。

5个降水个例24小时降水预报各个阈值的平均FSS评分

本章小结

从5个降水预报结果可以看出，在同化GPS ZTD资料后对降水强度和位置都一定程度的改进。在采用不同观测误差同化GPS ZTD资料对预报会产生较大影响。采用0.8cm为观测误差同化GPS ZTD资料能够有效改进降水位置预报，减小过报。采用1.99cm作为观测误差同化GPS ZTD资料其降水预报和控制试验相近，效果劣于0.8cm。所以采用0.8cm作为观测误差是更合适的选择。