一种局地集合资料同化方法
阅读说明：本技术 一种局地集合资料同化方法 (Method for assimilating local collection data ) 是由 王世璋 曾明剑 于 2021-09-27 设计创作，主要内容包括：本发明公开了一种局地集合资料同化方法,包括：选定分析时刻并生成分析时刻的模式背景场X～(f)和预报误差样本集合X'；读取模式背景场X～(f)、预报误差样本集合X',以及同化分析需要且能够获得的观测y～(o),并根据经纬度对模式预报区域和观测进行分区,并标记存储观测周围的模式格点的位置信息；根据计算条件和实际需求选定单点、单柱、或多柱组合的分析方式；对背景场X～(f)涉及的所有点、或单柱、或多柱组合进行同化分析；输出分析场。本发明的方法实现了模式变量背景误差协方差在同化过程中的直接计算,避免了直接使用集合样本而导致的分析精度下降的问题(次优解),同时解决了在同化垂直积分型观测时难以进行局地化的问题。(The invention discloses a method for assimilating locally-collected data, which comprises the following steps: mode background field X for selecting analysis time and generating analysis time f And a prediction error sample set X'; read mode background field X f A set of forecast error samples X', and observations y needed and available for assimilation analysis o Partitioning the mode forecasting area and observation according to the longitude and latitude, and marking and storing the position information of mode lattice points around the observation; selecting a single-point, single-column or multi-column combined analysis mode according to the calculation conditions and the actual requirements; to the background field X f All points involved, or single columns, or multi-column combinations, are subjected to assimilation analysis; and outputting an analysis field. The method of the invention realizes the direct calculation of the covariance of the background error of the mode variable in the assimilation process, and avoids the analysis precision caused by directly using the set sampleThe problem of degradation of degree (suboptimal solution) and the problem of difficulty in localization in assimilation vertical integration observation are solved.)
技术领域
本发明涉及天气预报技术，具体涉及一种局地集合资料同化方法。
背景技术
目前，现代天气预报的水平与数值天气预报的准确程度紧密相关，而数值天气预报的效果又依赖于初值的准确程度。当前，提高初值准确程度的主要方式是资料同化。欧洲中期天气预报中心（ECMWF）在数值天气预报的领先在相当程度上是其先进资料同化系统的贡献。
提高模式初值的准确程度是大气资料同化的目的，其途径是根据背景误差协方差，使用数学方法，在观测和模式初值之间进行拟合，使拟合得到的分析场在理论上更接近大气的真值。因此，资料同化自身的发展也主要围绕如何进行拟合，如何构建和使用背景误差协方差展开。
经过几十年的发展，资料同化的拟合方法，从早期的逐步订正、最优插值，发展到了到当前几种主流的方法：变分同化和集合卡尔曼滤波同化。资料同化拟合方法的出发点，可以是最小二乘，也可以是贝叶斯估计，但最终都可以表示为一个包含观测与背景信息的，以观测误差协方差和背景误差协方差为权重的代价函数。资料同化的分析场是使得这个函数达到极小值的解。资料同化的背景误差协方差，从早期的距离相关发展到当前包含多变量平衡的热动力统计相关。与早期相比，当前的主流资料同化方法能够直接同化高时空分辨率的非常规观测，能够考虑模式预报变量间的交叉误差协方差，其分析场具有更高的准确度（主要来自高分辨率非常规观测资料的贡献）并具有更高的热动力协调性（来自热动力平衡方程约束和流依赖背景误差协方差的贡献）。近年来，资料同化发展的热点之一是混合同化，混合同化是变分同化与集合卡尔曼滤波同化基于给定权重系数的结合，其特点在于兼顾了变分同化的稳定效果（静态背景误差协方差）与集合卡尔曼滤波流依赖背景误差协方差对快速变化误差的跟踪能力。
当前全球主要业务部门使用的资料同化方法基本都为变分同化、集合卡尔曼滤波同化以及混合同化，包括欧洲中期预报中心、美国环境气候预报中心、日本气象局、英国气象局、法国气象局、加拿大气象局、中国气象局等。其中欧洲中期预报中心、日本气象局、英国气象局、中国气象局等主要使用变分方法，加拿大气象局使用集合卡尔曼滤波方法。混合同化是上述主要气象业务部门都在发展的技术，尽管实现途径不同，但均考虑了静态背景误差协方差与流依赖背景误差协方差的结合。
资料同化发展的另一个驱动力与数值计算条件和数值预报模式、观测手段的发展密切相关。当计算条件不足时，资料同化方法需要进行很多假设和简化以满足计算条件限制，即使先进的同化方法也难以得到理想的结果；而计算条件充足时，不再需要假设和简化，同化效果也能够得到保证。当前，随着计算条件的不断提升，数值模式一方面为提高计算精度，使用越来越高的分辨率，网格数量也越来越多；另一方面，为适应不同的应用场景和需求，模式网格结构种类越发丰富，不再局限于正四边形网格，三角形、六边形网格，非结构变分辨率网格等多种结构在数值预报模式中得到了应用。因此，适应超大网格数量规模且不受网格结构限制成为未来资料同化方案发展的趋势。此外，大气观测观测手段不断增加，大量不同时空分辨率的观测资料也对同化方案提出了多源多尺度观测同时同化的需求。
在面对适应超大网格数量规模且不受网格结构限制的发展趋势，以及多源多尺度观测同时同化的需求时，现有同化方法存在一些不足，可以通过本发明方法得到较为理想的解决。
a、现有技术通用表达形式
现有同化方法都可以用一个代价函数J表示(参考Gao et al., 2004的版本)：
， （1）
其中v是控制变量，C是背景误差协方差的平方根，d是观测与投影到观测变量空间的模式初值的差值，H是上述投影的切线性算子，R是观测误差协方差矩阵。资料同化的目标就是求解使得J的梯度等于0 的控制变量v的解。当J的梯度等于0时，得到表达式
， （2）
其中I是单位矩阵。
变分方法、集合卡尔曼滤波方法以及混合方法的主要区别在于如何求解（2）式。
b、传统变分方法
在主流变分同化系统（如WRFDA(Barker et al., 2012)、GSI(Hu et al., 2017)、ARPS3Dvar(Gao et al., 2004)等）中，（2）式通过下降算法求解，例如共轭梯度法。首先需要根据统计信息重构背景误差协方差矩阵的平方根C。由于完整模式空间的背景误差协方差矩阵非常大（107×107），无法直接统计计算，因此通常会使用一些函数进行构造C矩阵。C通常表示为三个部分的乘积：Cp、Cv、Ch，其中Cp是模式变量与控制变量之间的物理变换，Cv是垂直误差协方差矩阵，Ch是水平误差矩阵。Cp因为往往涉及差分方程，所以与模式网格结构紧密相关。不同网格分辨率、不同网格结构都需要重新构建Cp。Cv通常通过特征值分解来实现计算。Ch在主流变分同化系统中均通过递归滤波方法实现计算。
上述三个C的分量中，与本发明的提出有密切联系的是Ch所对应的递归滤波方法，下面将予以详述。其余两个分量，因具体实现形式种类繁多，且与本发明关联度较弱，故不予赘述。根据递归滤波器的定义(Gao et al., 2004)：
, （3）
其中ψ i 是第i个格点的初始值（例如某步迭代中v在第i个点的值）, ϕ i 是从i到n的滤波之后的数值， χ i 是一次完整递归滤波后的数值，β决定了滤波的去相关系数：
, （4）
其中L 是去相关尺度, ∆s是任意一个方向的空间网格距 (例如, ∆s在东西方向上是∆x 而在南北方向上是∆y), 递归滤波通过次数为N pass。上述递归滤波器通常先在模式预报区域的东西方向自西向东依次计算，再在南北方向上从南往北依次计算。因此需要使用规则的四边形网格。
在确定了C的计算方案后，变分同化系统使用下降算法求解v。这里以共轭梯度法(Shewchuk 1994)为例。首先定义变量：
， （5）
然后以迭代的形式计算
。 （6）
当小于一定的阈值，如10-6时，停止计算，输出此时的vi+1值。模式变量空间的分析增量通过计算得到。
需要注意的是，变分方法通常是在整个模式预报区域上进行求解。
c、集合卡尔曼滤波方法
集合卡尔曼滤波方法有许多变种，其中具备高效并行计算能力和同时同化多种观测资料能力的方法是LETKF方法（Hunt.et al., 2007）。该方法与本发明关系也最为密切，此处予以详述。其他变种如EnSRF、EAKF、ETKF等在此不予赘述。LETKF方法通过集合的方式估计背景误差协方差，因此（2）式的C在LETKF中不再表示背景误差协方差矩阵的平方根，而是集合成员X与集合平均的差。
首先定义变量
。 （7）
然后对的转置进行SVD奇异向量分解，得到，其中E是左奇异向量矩阵，Σ是奇异值矩阵，D是右奇异向量矩阵。将SVD分解结果代入（2）式得到：
。 （8）
其中I是单位矩阵。模式变量空间的分析增量通过计算得到。
LETKF在串行版本中会依次更新每个格点上的模式变量，但每个格点的计算相互独立，因此可以高度并行化。此外，由于是单个点计算，因此其局地化方案可以通过调整观测误差协方差R来实现。相对于LETKF正在分析的格点，远距离的观测给予较大的误差值，近距离的观测给予相对较小的误差值。局地化调整系数根据距离相关函数确定。该函数在接近当前分析模式格点的区域数值接近1.0，在远离当前分析模式格点的区域数值趋于0.0，函数的分布形态接近高斯分布。
d、混合同化方法
现有混合同化方法的核心是静态背景误差协方差矩阵和流依赖背景误差协方差矩阵的结合（Wang et al., 2007）。有两种形式的混合：
， （9）
和 ， （10）
其中B表示背景误差协方差矩阵，下标hybrid、s、d分别表示混合背景误差协方差、静态（static）背景误差协方差和动态(dynamic)背景误差协方差，x表示模式变量空间的向量，是混合权重系数。在混合同化方面，本发明直接采用（9）式的形式，不做任何创新。
上述现有技术存在不足的一个主要原因在于背景误差协方差的构造。传统变分方法使用递归滤波构造水平误差协方差导致难以适应不同结构的网格，而LETKF在集合空间进行最小化，使得计算精度相对较低，且必须依赖集合，难以处理路径积分型观测。因此，如何构造背景误差协方差矩阵成为解决上述问题的关键。
发明内容
本发明的主要目的在于提供一种局地集合资料同化方法，较之变分同化方法，不受模式网格结构影响、能够同时同化多源多尺度观测资料；较之LETKF方法，能够更好的处理路径积分型观测，计算效率和计算精度更高，易于实现静态与动态背景误差协方差的组合。面向大规模计算需求，兼容多种网格类型，具备多源多尺度观测同时同化能力，并在一定程度上兼顾同化的计算精度。
本发明采用的技术方案是：一种局地集合资料同化方法，包括：选定分析时刻并生成分析时刻的模式背景场Xf和预报误差样本集合X'；读取模式背景场Xf、预报误差样本集合X'，以及同化分析需要且能够获得的观测yo，并根据经纬度对模式预报区域和观测进行分区，并标记存储观测周围的模式格点的位置信息；根据计算条件和实际需求选定单点、单柱、或多柱组合的分析方式；对背景场Xf涉及的所有点、或单柱、或多柱组合进行同化分析；输出分析场。
进一步地，所述根据经纬度对模式预报区域和观测进行分区包括：根据计算条件和需求，设定分区的块数，东西方向pnx块，南北方向pny块；遍历背景场第一层的每个格点，根据其经纬度标记其在分区中的位置；遍历每一个观测，根据其经纬度标记其在分区中的位置；遍历每一个观测，搜索每一个观测周围最近的n个模式格点，并记录其编号。
更进一步地，所述对背景场Xf涉及的所有点、或单柱、或多柱组合进行同化分析包括以下步骤：根据选定的观测yo确定背景误差协方差矩阵涉及的模式变量，并根据分区信息以及预先设定的局地化距离确定影响当前点、或单柱、或多柱组合的观测子集；将确定的模式变量及其扰动插值到确定的观测子集对应的位置，分别存储到HiXf和HiX'，Hi是模式网格点到观测网格点的插值算子；计算HiX’中每一个模式变量与Hi X'中其他模式变量的相关系数，得到各变量的标准偏差S和观测网格上的模式变量相关系数HiCcorrel，记为；计算各元素之间的局地化系数，得到局地化系数矩阵ρ，并与做元素-元素的乘积（Schur乘积），得到局地化后的；若中的元素涉及不同的模式变量，则对该元素乘以ρm，ρm的范围是0~1.0；计算调整系数λ，根据，其中trace（）表示计算矩阵的迹；计算；计算观测对应的innovation（d=yo-H oXf），其中H o是将模式变量投影到观测变量的观测算子，可以是非线性算子，H o的内容依具体观测而变；计算和，其中Ho是H o的切线算子；根据共轭梯度法求解方程得到控制变量v；计算选定的模式网格点上的变量xm的扰动Xm’与中各个变量的协方差以及局地化系数，以元素-元素的乘积得到局地化后的行向量Cm，Cm与v相乘，再乘以调整系数λ，得到xm的增量；计算得到xm的分析。
更进一步地，所述对背景场Xf涉及的所有点、或单柱、或多柱组合进行同化分析还包括：在循环代码之前，加入OpenMP引导语句，将关于每个点、或单柱、或多柱组合的编号设为私有变量，每个点、或单柱、或多柱组合的同化分析所用的观测数量、观测子集、误差协方差、控制变量等设为私有变量，其余设为共有变量；在循环代码结束后，加入OpenMP引导语句，结束OpenMP并行。
本发明的优点：本发明的局地集合资料同化方法。该方法构建了一个投影到观测网格的模式变量的局地背景误差协方差，实现了模式变量背景误差协方差在同化过程中的直接计算，避免了传统局地集合同化方法直接使用集合样本而导致的分析精度下降的问题（次优解），同时解决了传统局地集合同化方法在同化垂直积分型观测时难以进行局地化的问题。
所述局地集合资料同化方法，在相对传统局地集合同化方法提高计算精度的同时，引入了共轭梯度下降算法，以及单柱、多柱组合分析的方法，能够大幅减少观测重叠区域的重复计算，显著加快同化分析的速度，为该方法未来在业务或实际使用中提供了计算效率保障。
所述局地集合资料同化方法，相对传统局地集合同化方法，能够结合静态背景误差实现混合同化，有利于提高在预报误差样本数量不足时的同化分析效果。
所述局地集合资料同化方法，相对于使用递归滤波技术的传统变分框架的同化方法，能够不受模式预报网格结构的限制，能够适应更多应用场景，并且易于移植。
除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。
附图说明
构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。
图1是本发明的局地集合同化方法与LETKF的计算效率对比图；
图2是本发明的模式变量(a)u,(b)v,(c)θ,(d)q v的背景误差图；
图3是本发明的同化观测为观测类型2的模式变量(a)u,(b)v,(c)θ,(d)q v的背景误差图；
图4是本发明的同化观测为观测类型3的模式变量(a)u,(b)v,(c)θ,(d)q v的背景误差图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
图1 局地集合同化方法与LETKF的计算效率对比图（黑色实线为线性回归拟合的趋势线），纵轴为局地集合同化方法完成所有格点分析所需墙钟时间与LETKF所用时间的比值，横轴观测总数；
图2 模式变量(a)u,(b)v,(c)θ,(d)q v的背景误差图（交叉标记实线）、局地集合同化方法（LEDA）分析误差（实线）和LETKF分析误差（虚线）的垂直分布。纵轴为模式垂直层，横轴为误差大小。同化观测为观测类型1；
图3是本发明的同化观测为观测类型2的模式变量(a)u,(b)v,(c)θ,(d)q v的背景误差图；
图4是本发明的同化观测为观测类型3的模式变量(a)u,(b)v,(c)θ,(d)q v的背景误差图。
本发明的局地集合资料同化方法，具体包括如下步骤：
I）选定分析时刻并生成分析时刻的模式背景场Xf和预报误差样本集合X'；
II）读取模式背景场Xf、预报误差样本集合X'，以及同化分析需要且能够获得的观测yo，并根据经纬度对模式预报区域和观测进行分区，并标记存储观测周围的模式格点的位置信息；
III）根据计算条件和实际需求选定单点、单柱、或多柱组合的分析方式（多柱组合分析消耗更多的内存，但同化分析速度更快）；
IV）基于III）的选择，根据同化公式（1）、（2）、（5）、（6），并利用局地背景误差协方差公式，对背景场Xf涉及的所有点、或单柱、或多柱组合进行同化分析；
V）完成步骤IV）后输出分析场。
所述步骤II）中，模式预报区域和观测的分区，包括如下具体步骤：
（1）根据计算条件和需求，设定分区的块数，东西方向pnx块，南北方向pny块；
（2）遍历背景场第一层的每个格点，根据其经纬度标记其在分区中的位置，例如pnx和pny都是10，则背景场最中间的点，及其临近点的分区编号均为（5，5）；
（3）遍历每一个观测，根据其经纬度标记其在分区中的位置；
（4）基于（2）、（3）的结果，遍历每一个观测，搜索每一个观测周围最近的n个模式格点（n根据用户需要给定，通常为5~9），并记录其编号。
所述步骤IV）中，每个点、或单柱、或多柱组合的同化分析，包括如下具体步骤：
（1）根据II）中选定的观测yo确定中矩阵涉及的模式变量（包括涉及的垂直层数）；并根据II）中的分区信息以及预先设定的局地化距离确定影响当前点、或单柱、或多柱组合的观测子集；
（2）将（1）中确定的模式变量（背景场或集合平均）及其扰动（集合成员与集合平均的差值）插值到（1）中所确定的观测子集对应的位置，分别存储到HiXf和HiX'（Hi是模式网格点到观测网格点的插值算子，插值算法由用户确定，本发明以距离权重插值为例）；插值使用的模式变量的位置由步骤II）确定；
（3）计算HiX’中每一个模式变量（即每一行对应的模式变量）与Hi X'中其他模式变量（其他行）的相关系数，得到各变量的标准偏差S和HiCcorrel（记为）；
（4）计算各元素之间的局地化系数，得到局地化系数矩阵ρ，并与做元素-元素的乘积（Schur乘积），得到局地化后的；若中的元素涉及不同的模式变量，则对该元素乘以ρm，ρm的范围是0~1.0；优选的，ρm=0.5；
（5）计算调整系数λ，根据，其中trace（）表示计算矩阵的迹；
（6）计算
（7）计算观测对应的innovation（d=yo-H oXf），其中H o是将模式变量投影到观测变量的观测算子，可以是非线性算子，H o的内容依具体观测而变；
（8）计算和，其中Ho是H o的切线算子；
（9）根据共轭梯度法的（5）式和（6）式计算得到控制变量v；
（10）计算选定的模式网格点上的变量xm的扰动Xm’与中各个变量的协方差以及局地化系数，以元素-元素的乘积（Schur乘积）得到局地化后的行向量Cm，Cm与v进行点乘，再乘以调整系数λ，得到xm的增量；
（11）计算得到xm的分析。
所述步骤IV）中，对所有点、或单柱、或多柱组合进行同化分析，使用如OpenMP等方法进行并行计算，提高计算效率，包括如下步骤：
（1）在步骤IV）的循环代码之前，加入OpenMP引导语句，将步骤IV）中关于每个点、或单柱、或多柱组合的编号为私有变量，每个点、或单柱、或多柱组合的同化分析所用的观测数量、观测子集、误差协方差、控制变量等设为私有变量，其余设为共有变量；
（2）正常执行步骤IV）；
（3）在步骤IV）循环代码结束后，加入OpenMP引导语句，结束OpenMP并行。
以下给出一段具体实施例，该例不需要进行本地数值预报，使用模拟观测，便于用户在计算条件有限且没有实际观测的情况下对同化算法测试检验：
(I.1）从NCDC官网https://www.ncdc.noaa.gov/data-access/model-data下载分辨率为0.5°的GFS数据，其中下载任意一天的12UTC的36小时预报作为模式背景场Xf。下载该时刻前5天的，预报到该时刻的GFS预报数据用于生成预报误差样本集合X'（时间滞后法生成集合成员）。下载该时刻的GFS分析场作为真值（用于检验分析结果）。使用WRF v3.9.1模式和WPS v4.1将下载的GFS数据转成区域数据，网格分辨率、区域范围可由用户自行设定。这里以15km网格分辨率和150×150网格点数为例。
(I.2)生成模拟观测：
模拟观测类型1：常规探空观测，假定Ho=H o=I，其中I是单位矩阵。观测变量为水平风场的两个分量u和v，位温θ，水汽混合比q v、以及地面气压ps。类型1模拟观测的东西和南北方向密度均为每7个格点一个，观测误差u和v为0.5 m/s，位温θ为0.5 K，水汽混合比q v为0.0005 kg/kg、以及地面气压ps为10 Pa。除了ps，其余观测变量在垂直方向上均为每2模式层取一个。
模拟观测类型2：整层可降水水汽含量（pwv）。Ho=H o=1/gΣ（qv*dp），Σ表示从模式第一层到最高层的积分，dp表示每层之间的气压差，qv是每层的水汽混合比。观测密度为东西和南北方向每4个格点一个，观测误差数值大小1.0。
模拟观测类型3：风廓线观测，包含两个变量：风向和风速。风向的H o=arctan(u/v)+nπ，其中n是整数，使得H o在正负0.5π之间；风向Ho=u/sδv-v/sδu。风速的H o=（u2+v2）0.5；风向Ho=u/sδu+v/sδv，其中s=（u2+v2）0.5。观测密度为东西和南北方向每3个格点一个，垂直方向包括所有模式层，风向和风速观测误差大小分别为0.5 m/s和0.1弧度。
（I.3）设定同化系统运行参数，包括三种类型观测的局地化半径和搜索距离等。
（I.4）读取模式背景场Xf和预报误差样本集合X',和(I.2)所述三种观测。
（I.5）进行模式数据和观测的经纬度分区，并标记观测周围的模式格点：
将预报区域的东西和南北方向各分出10块，即pnx和pny等于10；遍历背景场第一层的每个格点，根据其经纬度标记其在分区中的位置；遍历每一个观测，根据其经纬度标记其在分区中的位置；遍历每一个观测，搜索每一个观测周围最近的5个模式格点，并记录其编号，搜索仅在相邻分区中进行。
（I.6）若使用多柱组合分析方式，则记录分组数量和每组对应的模式网格点（本例中使用9柱组合分析）。
（I.7）对每一个多柱组合进行同化分析：
根据II）中的分区信息以及预先设定的局地化距离确定影响当前多柱组合的观测子集；将模式变量（背景场或集合平均）及其扰动（集合成员与集合平均的差）插值到观测子集中观测网格的位置，分别存储到HiXf和HiX'（Hi是模式网格点到观测网格点的距离权重插值算子，其中参与插值的模式变量位置由（I.5）确定；计算HiX’中每一个模式变量（即每一行对应的模式变量）与Hi X'中其他模式变量（其他行）的相关系数，得到各变量的标准偏差S和HiCcorrel（记为）；计算各元素之间的局地化系数，得到局地化系数矩阵ρ，并
与做元素-元素的乘积（Schur乘积），得到局地化后的；若中的元素涉及不同的模式变量，则对该元素乘以ρm（=0.5）；计算调整系数λ，根据；计算计算观测对应的innovation（d=yo-H oXf），其中H o是将模式变量投影到观测变量的观测算子，可以是非线性算子，本例中H o的内容依具(I.2)确定；计算和，其中Ho是H o的切线算子，由(I.2)确定；根据共轭梯度法的（5）式和（6）式计算得到控制变量v；计算当前多柱组合包含的每一个模式网格点上的变量xm的扰动Xm’与中各个变量的协方差以及局地化系数（包括ρm），以元素-元素的乘积（Schur乘积）得到局地化后的行向量Cm，Cm与v进行点乘，再乘以调整系数λ，得到xm的增量；计算得到xm的分析；进入下一个多柱组合（若开启OpenMP并行模式，则根据使用的CPU核数确定同时计算的多柱组合数量n，在I.7中同时更新n个多柱组合）。
（I.8）完成所有模式网格点和模式变量的同化分析，输出分析结果，结束同化运行流程。
本发明所设计的局地集合资料同化方案解决的关键技术问题细节如下：
1、通过空间投影，实现局地背景误差协方差矩阵（平方根）的直接计算。
为解决背景误差协方差矩阵难以直接计算的问题，本发明对（2）式中的背景误差协方差的平方根C矩阵进行了重新的设计。通常（2）式中的C矩阵涉及完整预报网格区域的所有模式变量，因而自由度非常高。而本发明对C矩阵进行了投影。（2）式中的HC可以分解为两个部分Ho和HiC，其中Ho和Hi分别为模式变量向观测变量的转换与模式网格向观测网格的投影（插值）。当观测涉及整层积分时（如pwv观测），Hi包含多个气柱的整层投影（插值）。这一处理方式解决了LETKF难以处理路径积分型观测的问题。本发明使用HiC作为背景误差协方差矩阵的平方根，记为。与C比较可以发现，的元素数量大幅减少，因为只涉及计算Ho所需要的模式变量。同时，参考LETKF进行局地化的思路：远离当前分析格点的观测可以不予考虑，则对于任意一个分析格点而言，只涉及到一个局地范围内的模式变量，因此可以直接显式计算。这里将称为局地背景误差协方差。当使用集合方式估计C时，理论上可以先计算局地的B=XXT，并通过Schur product将局地化系数矩阵作用到B上，再通过SVD分解得到经过局地化的。但是该方法计算量较大，计算效率较低，因此本发明设计了一个直接估计矩阵的方法：
， （11）
其中ρ是局地化系数矩阵，Ccorrel是背景误差的相关系数矩阵，S是一个对角阵，元素是模式变量的标准偏差，λ是一个量级调整系数，使得的迹与Ccorrel的迹相等（即两个矩阵的量级相当）。若不乘以这个调整系数，则（11）式估计的背景误差协方差矩阵的特征值为原背景误差协方差矩阵B（B=SCcorrelST）特征值的平方倍，背景误差协方差过大，使同化分析结果不合理。λ的计算表达式为,其中N是Ccorrel的秩（或行数），trace（）表示计算矩阵的迹。
这里对矩阵的大小做一个说明。假定Ho是一个nobs×nh的矩阵，nobs是一个模式格点周围一定半径范围内的观测数量，nh是观测涉及的模式变量数（格点数与模式变量种类数的乘积），则的大小是nh×nh。实际计算中，nh的量级大约在100~103，取决于局地观测的数量和观测涉及的模式变量数量，远小于完整C矩阵的107量级，完全在当前单核计算能力范围内。由于矩阵经过了局地化处理，因此是一个以对角元素为主的矩阵，其条件数（condition number，最大特征值与最小特征值的比值）远小于LETKF的X，因此病态程度较低。而且，由于投影过程并未进行任何假设，因此背景误差协方差的秩并没有减少。
此外，由于的元素可以直接通过集合样本统计得到，局地化也可以直接根据观测间的经纬度差异计算，因此使用的局地集合资料同化方法可以避免如递归滤波方法受网格结构的限制。
直接计算背景误差协方差的另一个优势是该方法为调节模式变量间协方差大小的提供了可能。由于数值模式具有高度非线性，模式变量间的协方差往往难以准确表征两个模式变量之间的变化趋势。若协方差正负有误，或者量级过大，则会产生错误或过大的增量，导致同化分析不能达到降低误差的效果。适度调节变量间协方差的大小，有利于减少过大且不合理的分析增量，有利于提高同化分析效果。具体方式是：记录Ccorrel中每个元素对应的两个模式变量的类型，若两个变量类型不一样，则乘以一个预先设定的系数ρm。ρm的大小在0~1之间。在本发明中，建议使用0.5。
2、使用共轭梯度法和多气柱同时分析的方法实现同化分析的高效快速计算。
共轭梯度法所用（5）式和（6）式需要矩阵对称且正定。实测当使用LETKF的方式（集合样本）估计C时，的条件数非常高，共轭梯度法求解无法收敛。因此只能通过SVD方式就行求解。而本发明所用矩阵（11式）由于条件数大幅降低，因此构成的矩阵能够使共轭梯度法的求解过程收敛。共轭梯度法由于本身计算以向量的点乘为主，因此计算速度远快于SVD分解的速度，尤其是面对大背景误差协方差矩阵的情况。
由于相当于一个在水平方向上稀疏化的背景误差协方差矩阵的平方根，因此其局地化过程并不针对某个模式格点，因此可以选择单点、单柱（一个气柱）、多柱（多个气柱）进行同时分析，减少了观测重叠区域的重复计算，减少了观测重叠区域的重复计算，加快了同化分析的速度和运行效率。
3、使同化方案易于实现不同类型协方差的混合
与LETKF方法无法使用静态背景误差协方差相比，本发明设计的局地集合资料同化方法可以很方便的根据（9）式实现静态与动态背景误差协方差的混合。当使用的静态背景误差协方差以距离相关函数估计时，的元素可以表示为：
， （12）
其中下标i,j表示矩阵的第i行第j列，上标Ens表示集合样本，ρ表示距离相关系数，S表示第i行对应的模式变量的标准偏差。当权重系数γ减为0时，变为完全的静态协方差。
实施例
（1）相对使用矩阵分解（特征值分解或SVD分解）的LETKF，本发明具有更高的计算速度
这里给出一组测试对比，测试考虑了不同的观测密度。从图1可以看出本发明所述局地集合同化方法在计算用时上远小于LETKF方法，同化大约15万观测需要的时间仅为LETKF的20%或更低。产生这一结果的原因，一方面是由于本发明所述局地集合同化方法减少了重复计算并应用了共轭梯度法，另一方面是多气柱组合（如这里的9气柱组合）减少了同化的单元数量。在LETKF中，单元数量为格点总数；在本发明所述局地集合同化方法中，单元数量为格点总数除以多气柱组合包含的格点数。因此每个单元内的气柱越多，总单元数量就越少。但是每个单元涉及的观测也会增加，因此每单元内的气柱数量不能无线增加。使用9个气柱左右即可显著提升计算速度。
（2）相对LETKF，本发明所述局地集合同化方法的分析误差相当或较小
这里给出一组测试对比，测试考虑了不同的观测类型及组合（参考D2.2具体实施方式（I.2）生成模拟观测）。从图2中可以看出，在同化常规探空观测（类型1）时，本发明所述局地集合同化方法分析的u和v在低层的误差小于LETKF的误差。分析的θ在中低层误差均小于LETKF的误差。两种方法在水汽q v的分析上大致相当。从图3中可以看出，当同化观测类型2（可降水水汽含量pwv）时，本发明所述局地集合同化方法与LETKF对u和v的分析效果大致相当，但没有像LETKF一样造成低层θ误差的增加。两者对水汽的调整水平基本相当。从图4中可以看出，当同化非线性观测（风向、风速——观测类型3）时，两者的分析效果均不如同化直接观测（类型1），但是本发明所述局地集合同化方法的风场分量u和v在低层和高层略小于LETKF，仅在中层略大于后者。同时本发明所述局地集合同化方法避免了LETKF导致的非观测变量（θ和qv）的误差增加。整体而言，本发明所述局地集合同化方法在分析效果上至少有LETKF相当，对直接观测有更好的分析能力，对非直接观测的负影响（误差增大）较小。
（3）本发明所述局地集合同化方法具有较高的可扩展性和易移植性
与使用递归滤波技术的变分同化方法相比，本发明所述局地集合同化方法不受网格结构、分辨率连续变化等情况的限制，可应用于不同类型的数值模式；与LETKF方法相比，本发明所述局地集合同化方法能够使用混合背景误差协方差，甚至完全静态的背景误差协方差，可在集合数量有限甚至没有集合成员的情况下完成同化分析。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
- 上一篇：一种医用注射器针头装配设备
- 下一篇：水生植被的类型确定方法及装置