基于L1范数正则项约束的不连续资料三维/四维变分融合研究

引用本文

王根, 盛绍学, 刘惠兰, 吴蓉, 杨寅. 基于L1范数正则项约束的不连续资料三维/四维变分融合研究. 地球科学进展, 2017, 32(7): 757-768
Wang Gen, Sheng Shaoxue, Liu Huilan, Wu Rong, Yang Yin. Discontinuous Data 3D/4D Variation Fusion Based on the Constraint of L1 Norm Regularization Term. Advances in Earth Science, 2017, 32(7): 757-768 复制到剪切板

Doi:10.11867/j.issn.1001-8166.2017.07.0757
Permissions

基于L1范数正则项约束的不连续资料三维/四维变分融合研究

王根^1,², 盛绍学¹, 刘惠兰¹, 吴蓉³, 杨寅⁴

1.安徽省气象信息中心安徽省大气科学与卫星遥感重点实验室, 安徽合肥 230031

2.中国气象局沈阳大气环境研究所, 辽宁沈阳 110000

3.安徽省气候中心, 安徽合肥 230031

4.国家气象中心, 北京 100081

作者简介:王根(1983-),男,江苏泰州人,工程师,主要从事卫星资料同化、GRAPES数值模拟和多源数据融合研究.E-mail:203wanggen@163.com

基金:安徽省自然科学基金项目“广义变分同化AIRS水汽通道亮温及在安徽强对流天气预报中的应用研究”(编号:1708085QD89); 淮河流域气象开放研究基金项目“基于地面和卫星观测反演的江淮流域降水资料融合算法研究”(编号:HRM201407)资助

摘要

经典三维/四维变分融合基于误差服从高斯分布,在极小化迭代时涉及到求解目标泛函梯度,若资料不连续则不可微,从而无法求解相应梯度,故理论要求所融合的资料必须具有“连续性”。采用扩展经典三维/四维变分融合方法,显式地基于L1范数把先验知识作为正则项约束项耦合到经典变分融合模型。在实施过程中把资料映射到小波域,采用新的融合模型在“小波空间”完成资料融合后,再采用小波逆变换映射回“观测空间”。通过线性平流扩散方程作为四维预报模式进行理想试验,试验设计融合背景和观测资料不连续,即在某些点左右导数不相等,试验结果表明文中采用的方法可行。进一步将该方法用于多源降水资料融合试验,采用基于GAMMA拟合函数的概率密度匹配法(Probability Density Function matching method,PDF)进行CMORPH反演降水资料订正,再将订正后的资料与地面站观测资料进行融合。通过与参考场结构相似性度量,得到该方法能更好地保留代表一些天气现象的“离群点”。该融合方法为不连续资料融合,尤其是“跳变点”的变分融合奠定了理论基础并提供了可借鉴的方法。

关键词: 不连续资料; 变分融合; L1范数; 正则项; 小波空间

中图分类号:P468 文献标志码:A 文章编号:1001-8166(2017)07-0757-12

Discontinuous Data 3D/4D Variation Fusion Based on the Constraint of L1 Norm Regularization Term

Wang Gen^1,², Sheng Shaoxue¹, Liu Huilan¹, Wu Rong³, Yang Yin⁴

1.Anhui Meteorological Information Centre Anhui Key Laboratory of Atmospheric Science and Satellite Remote Sensing, Hefei 230031, China

2.The Institute of Atmospheric Environment, China Meteorological Administration, Shenyang 110000, China

3.Anhui Climate Center, Hefei 230031, China

4.National Meteorological Center of China, Beijing 100081, China

First author:Wang Gen(1983-),male,Taizhou City, Jiangsu Province, Engineer. Research areas include satellite data assimilation, numerical simulation of GRAPES and multi-source data fusion.E-mail:203wanggen@163.com

Fund:Project supported by the Natural Science Foundation of Anhui Province“Generalised variational assimilation of AIRS water vapor channel brightness temperature and the application study in severe convective weather in Anhui Province”(No.1708085QD89); The Open Research Fund of Huai River Basin in Meteorological “Study of precipitation data fusion algorithm in Jianghuai Basin based on the ground and satellite observations”(No.HRM201407)

Abstract

Classical 3D/4D variation fusion is based on the theory that error follows Gaussian distribution. When using minimization iteration, the gradient of objective function is involved, and the solution of which requires the continuity of data. This paper adopted the extended classical 3D/4D variation fusion method, and explicitly applied the prior knowledge, which was based on L1-norm, as regularization constraint to the classical variation fusion method. Original data was firstly projected into the wavelet domain during the implementation process, and new fusion model was adopted for data fusion in wavelet space, then inverse wavelet transform was used to project the result to the observation space. Ideal experiment was carried out by using linear advection-diffusion equation as four-dimensional prediction model, which made a hypothesis of the discontinuity with the data between background and observation, and that meant the derivatives between left and right were not equal on some points. The result of the experiment showed that the method adopted here was practicable. A further research was also done for multi-source precipitation fusion. Firstly, CMORPH inversion precipitation data were corrected through PDF (Probability Density Function, PDF) matching method based on GAMMA fitting function. Then corrected data was fused with the observation one. By comparison with the reference field, the result showed that this method can keep some outliers better, which might represent certain weather phenomenon. The L1-norm regularization variation fusion in this paper provided a possible way to deal with discrete data, especially for jump point.

Keyword: Discrete data; Variation fusion; L1-norm; Regularization term; Wavelet space.

Show Figures

1 引言

多源融合技术在提高水文— 气象预测能力、气候和海洋系统建模等方面发挥了核心作用^[1]。其基本思想是不断融入状态变量(如风、温度和压力等)观测值到物理状态中(如云和降水等微物理过程), 随时间进行迭代通过变分融合系统减少状态变量的不确定性^[2]。数据融合方法通常利用观测结果更新当前先验模型状态估计(背景场)进而产生统计意义上的“ 最优” 后验状态(分析场), 并用于预测下一个时间的物理状态^[3]。目前已有数值天气预报中心已经或正在建立变分融合系统。虽然各中心在建立时设计方案各有不同, 但原理上融合方案都要最小化一个平方型的评价函数, 即最小平方估计^[4]。最小平方估计对远离值非常敏感, 即使仅有几个边远值, 也会使最小平方估计结果完全脱离实际。文中把目前国际上通用融合方法记为“ 经典变分融合” , 常用方法分为两大类:一是群递归滤波方法, 其本质是把随时间演化得到的统计特性融合到系统以追踪最优状态(如基于卡尔曼滤波驱动的数据融合方法)^[5~7]; 二是变分法, 如三维/四维变分(3DVar/4DVar), 其本质是依靠模式直接优化每一瞬间所有可用的观测资料^{[8, 9]}。经典变分资料融合通常需要解决“ 光滑” 优化问题, 其解决方案是通过定义一个最小加权欧氏距离, 即观测和背景与真值的差值估计^[10], 权重由背景和观测误差协方差矩阵给定。从估计理论出发, 该问题等价于高斯噪声状态下的最大似然(Maximum Likelihood, ML)估计^[11]。经典3DVar/4DVar变分融合基于误差服从高斯分布, 在极小化迭代涉及到目标泛函梯度求解, 若资料不连续则不可导, 从而无法求解相应梯度, 故理论要求所融合的资料必须具有“ 连续性” ^[12]。

但实际上很多资料具有“ 不连续” 性, 从而引入了一个新的问题, 即如何进行“ 不连续” 资料融合?为此, 本文在Ebtehaj等^[13]研究的基础上, 显式地把分析状态底层结构信息, 即初始背景场作为正则项耦合到经典变分融合框架中。并在理论分析算法可行的基础上, 开展理想试验和不连续性多源降水资料融合试验。目前国内外关于多源降水资料融合主要分为“ 二源” ^[14~16]和“ 三源” 融合^[17]。“ 二源” 融合研究较多的是融合CMORPH反演降水资料与地面站观测资料, 采用概率密度匹配法(Probability Density Function matching method, PDF)和最优插值(Optimal Interpolation, OI)^[14]。因OI的数理理论要求资料具有连续性, 且OI是基于全局二次曲面最优, 易出现过度光滑一些“ 离群” 但正确的值, 从而导致结果易出现负值, 需进行后处理。“ 三源” 融合, 即融合地面— 雷达— 卫星反演资料。其具体思路为先采用PDF方法订正雷达和卫星估测降水的系统误差, 再采用贝叶斯融合方法将三者进行融合^[17]。

本文基于Ebtehaj等^[13], Tikhonov等^[18]和王根等^[19]的研究方法得到新的变分融合法, 为“ 不连续” , 如降水、云和大气可降水量(Precipitable Water Vapor, PWV)等一些“ 快变” 量融合奠定基础。该方法显式地包含“ 不连续” 资料底层分析状态的先验知识作为L1范数正则项^[13], 通过线性平流扩散方程进行理想试验(同化的背景资料和观测资料不连续, 且具有强的“ 突变点” )。并进一步将该方法用于经过PDF订正后的CMORPH反演降水资料与地面站观测资料融合, 以验证该方法是否能够提高分析状态的稳定性。引入正则化的方法为“ 连续” , 如气温、海洋表面温度( Sea Surface Temperature, SST)^[20]等一些“ 慢变” 量融合提供了一定借鉴。

2 经典3DVar/4DVar变分融合理论分析

经典3DVar/4DVar变分融合方法的核心为最小二乘估计理论^[4]。区别于OI和卡尔曼滤波等方法^[6], 变分具有灵活性, 主要体现在2个方面:①可以同时加入不同来源的观测资料; ②可以对目标代价函数加入约束项如正则项, 约束项可由拉格朗日乘子法加入到代价函数进行极小化迭代求得“ 数值解” 或者直接求“ 解析解” 。

假设在初始时间t₀未知真实状态是离散空间中每个包含m个元素的列向量x₀, x₀=[x_{0, 1}, …, x_0,_m]^T∈ R^m, 在时间间隔[t₀, …, t_k]观测和随机误差分别为y_i∈ Rⁿ和v_i, i=1, …, k, 其中n< < m。则观测与状态的关系为公式(1)^[13]:

y_i=H(x_i)+v_i (1)

式中:H:R^m→ Rⁿ表示非线性观测算子将状态空间映射到观测空间; v_i~N(0, R_i)是零均值和方差为R_i的高斯观测误差。变分融合的关键是针对不同观测资料构建不同的观测算子。对于多源降水资料融合, 观测算子H(· )可以假设为线性算子^[14]; 而对于卫星辐射资料, 该算子为辐射传输模式^{[19, 21]}。

3DVar/4DVar变分融合归结为在初始时刻极小化以下目标泛函^{[13, 22]}:

J₃_D_{, 4}_D(x₀, x₁, …, x_k)= $\begin{matrix} \overset{k}{\sum_{i = 0}} \end{matrix} \begin{matrix} (\frac{1}{2} ‖ y_{i} - H (x_{i}) ‖_{R_{i}^{- 1}}^{2}) \end{matrix}$ + $\begin{matrix} \frac{1}{2} \end{matrix}$ ‖ $\begin{matrix} x_{0}^{b} \end{matrix}$ -x₀ $\begin{matrix} ‖_{B^{- 1}}^{2} \end{matrix}$ (2)

x_i=M_0,_i(x₀), i=0, …, k (3)

式中:y_i为观测值, y_i∈ Rⁿ, i=0, …, k; $\begin{matrix} x_{0}^{b} \end{matrix}$ ∈ R^m和B∈ R^{m× m}分别表示背景状态和背景误差协方差矩阵。 $\begin{matrix} ‖ x ‖_{A}^{2} \end{matrix}$ =x^TAx标记为二次范数, 要求A是正定矩阵; 函数M_0,_i:R^m→ R^m标记为一个从初始状态t₀到t_i的非线性演变模型。当M=I时, 此时为3DVar, I为单位矩阵; 否则为4DVar。在数值天气预报中M为预报模式, 由一组复杂的预报方程(偏微分方程)构成。经典变分融合的目标是获得所谓最佳真实状态估计, 即分析状态 $\begin{matrix} x_{0}^{a} \end{matrix}$ ∈ R^m, 从变分融合角度看, 相当于最小化2个二次代价函数的总和。

定义M_0,_i的雅可比矩阵为M_0,_i。通过求解演变模型线性化后的切线性模型的伴随模式, 即可得到相应的雅克比矩阵^[19]。文中只考虑线性观测算子, 即H(x_i)=Hx_i, 则4DVar代价函数可以简化为下式:

J₄_D(x₀)= $\begin{matrix} \overset{k}{\sum_{i = 0}} \end{matrix} \begin{matrix} (\frac{1}{2} ‖ y_{i} - H M_{0, i} x_{0} ‖_{R_{i}^{- 1}}^{2}) \end{matrix}$ + $\begin{matrix} \frac{1}{2} \end{matrix} \begin{matrix} ‖ x_{0}^{b} - x_{0} ‖_{B^{- 1}}^{2} \end{matrix}$ (4)

假设 $\begin{matrix} \bar{y} \end{matrix}$ =[ $\begin{matrix} y_{0}^{T} \end{matrix}$ , …, $\begin{matrix} y_{k}^{T} \end{matrix}$ ]^T∈ R^N; N=n(k+1);

$\begin{matrix} \bar{H} \end{matrix}$ = $\begin{matrix} {[{(H M_{0, 0})}^{T}, \dots, (H M_{0, k})^{T}]}^{T} \end{matrix}$ ;

$\begin{matrix} \bar{R} \end{matrix}$ = $\begin{matrix} [\begin{matrix} R_{0} & 0 & \dots & 0 \\ 0 & R_{1} & \dots & ⋮ \\ ⋮ & ⋮ & 0 \\ 0 & \dots & 0 & R_{k} \end{matrix}] \end{matrix}$

则4DVar代价函数(4)式可进一步简化为极小化代价函数:

J₄_D(x₀)

= $\begin{matrix} \frac{1}{2} \end{matrix}$ ‖ $\begin{matrix} \bar{y} \end{matrix}$ - $\begin{matrix} \bar{H} \end{matrix}$ x₀ $\begin{matrix} ‖_{{\bar{R}}^{- 1}}^{2} \end{matrix}$ + $\begin{matrix} \frac{1}{2} \end{matrix}$ ‖ $\begin{matrix} x_{0}^{b} \end{matrix}$ -x₀ $\begin{matrix} ‖_{B^{- 1}}^{2} \end{matrix}$ (5)

公式(5)是状态x₀的二次光滑函数, 令导数为零, 则可得到分析状态:

$\begin{matrix} x_{0}^{a} \end{matrix}$ =( $\begin{matrix} {\bar{H}}^{T} \end{matrix} \begin{matrix} {\bar{R}}^{- 1} \end{matrix} \begin{matrix} \bar{H} \end{matrix}$ +B^-1)^-1( $\begin{matrix} {\bar{H}}^{T} \end{matrix} \begin{matrix} {\bar{R}}^{- 1} \end{matrix} \begin{matrix} \bar{y} \end{matrix}$ +B^-1 $\begin{matrix} x_{0}^{b} \end{matrix}$ ) (6)

3 基于小波域的资料“ 不连续” 变分融合理论分析

结合先验知识进一步约束“ 状态变量” 的变分融合, 得到改进的后验分析估计。根据最大后验(Maximum A Posteriori, MAP)估计和贝叶斯定理, 在4DVar代价函数中加入正则化约束项, 则有^{[13, 23]}:

$\begin{matrix} x_{0}^{a} \end{matrix}$ = $\begin{matrix} \underset{x_{0}}{argmin} \end{matrix}$ {J_R₄_D(x₀)+λ τ (x₀)} (7)

式中:τ (x₀)= $\begin{matrix} \sum_{i}^{k} \end{matrix}$ (Φ _ix_0,_i)^p=‖ Φ x₀ $\begin{matrix} ‖_{p}^{p} \end{matrix}$ 。本文将资料的不连续性引入到4DVar模型中, Φ ∈ R^{m× m}是“ 小波变换” , Lp范数定义为

‖ x‖ _p=(∑ |x_i|^p)¹^/p(p> 0)。

在公式(7)目标代价函数中, J和τ 分别是光滑和非光滑凸函数。在经典3DVar/4DVar变分融合时可以忽略τ 的属性, 但作为额外的先验知识可以用于提高融合结果的准确性。目前国际上用于解决大规模非光滑凸函数的优化问题主要有内点算法(interior point algorithms)^[24]和近端加速梯度方法(accelerated proximal gradient methods)^[25]等。本文中基于梯度投影算法(Gradient Projection, GP)进行迭代求解^{[13, 26]}。进一步得到:

$\begin{matrix} \{\begin{matrix} x_{0}^{a} = \underset{x_{0}}{argmin} {J_{R 4 D} (x_{0})} \\ ‖Φ x_{0} ‖_{p}^{p} \leq c \end{matrix} (8) \end{matrix}$

式中:c> 0, 通过约束Lp范数能得到更稳定的解。基于拉格朗日乘子法公式(8)的约束问题转化为求解无约束问题, 则有^[13]:

$\begin{matrix} x_{0}^{a} \end{matrix}$ = $\begin{matrix} \underset{x_{0}}{argmin} \end{matrix} \begin{matrix} \{\frac{1}{2} ‖ \bar{y} - \bar{H} x_{0} ‖_{{\bar{R}}^{- 1}}^{2} \end{matrix}$ + $\begin{matrix} \frac{1}{2} \end{matrix}$ ‖ $\begin{matrix} x_{0}^{b} \end{matrix}$ -x₀ $\begin{matrix} ‖_{B^{- 1}}^{2} \end{matrix}$ +λ ‖ Φ x₀ $\begin{matrix} ‖_{p}^{p} \end{matrix} \begin{matrix} \frac{}{}\} \end{matrix}$ (9)

λ 为拉格朗日乘子, 也称正则化参数, 公式(9)中要求λ 非负^[27]。当λ =0时, 公式(9)变成了经典变分融合法。较小的λ 值能够减小先验知识对结果的影响; 较大的λ 值能够引入更多的先验知识。因此λ 起着重要的平衡作用, 能同时保持足够地接近观测和背景状态。通常λ 参数的确定是通过统计交叉验证获得^[27]。

从概率角度出发, 公式(9)可以被视为正则化最大后验贝叶斯估计值。正则化约束是指状态量关于先验知识的概率分布^[28], 即p(x)∝ exp(-λ ‖ Φ x $\begin{matrix} ‖_{p}^{p} \end{matrix}$ )。

前期Wang等^[4]把稳健性较强的M— 估计(L2— 估计、Huber— 估计、Cauchy— 估计和Fair— 估计)耦合到经典变分方法中, 本文则把先验知识作为正则项。由于资料可能存在一些正确但可用的“ 离群点” , 为保留这些点, 本文基于前期Wang等^[4]的研究基础结合Ebtehaj等^[13]提出的方法, 将较稳健的L1范数耦合到4DVar作为目标泛函的正则项。并代替公式(9)中的Lp, 则有:

$\begin{matrix} x_{0}^{a} \end{matrix}$ = $\begin{matrix} \underset{x_{0}}{argmin} \end{matrix} \begin{matrix} \{\frac{1}{2} ‖ \bar{y} - \bar{H} x_{0} ‖_{{\bar{R}}^{- 1}}^{2} \end{matrix}$ + $\begin{matrix} \frac{1}{2} \end{matrix}$ ‖ $\begin{matrix} x_{0}^{b} \end{matrix}$ -x₀ $\begin{matrix} ‖_{B^{- 1}}^{2} \end{matrix}$ +λ ‖ Φ x₀‖ ₁ $\begin{matrix} \frac{}{}\} \end{matrix}$ (10)

公式(10)记为“ R4D” 。由L1代价函数性质决定其较适合融合一些“ 快变” 的变量如降水等。L1范数正则化4DVar目标泛函进一步定义为:

$\begin{matrix} \{\begin{matrix} x_{0}^{a} = \underset{x_{0}}{argmin} {‖ Φ x_{0} ‖_{1}} \\ J_{R 4 D} (x_{0}) \leq c \end{matrix} (11) \end{matrix}$

式中:c> 0, 此为求解带约束项的范数极小化问题。如果Φ 作为一阶偏导算子, 使用L2范数正则化(λ ‖ Φ x₀ $\begin{matrix} ‖_{2}^{2} \end{matrix}$ ), 得到的分析增量能量最小, 但会减少分析的可变性和对解的“ 过度” 平滑。平滑的L2范数正则化导数空间, 较适合连续和充分光滑的“ 状态变量” , 如气温。对分段光滑的状态变量或孤立奇异的“ 突变点” , 在导数空间较适合使用L1范数正则化(λ ‖ Φ x₀‖ ₁), 其能隐式地约束目标泛函并能防止在边缘和跳跃间断点产生额外平滑, 即减弱了经典变分融合的“ 过度” 拟合。

4 基于线性平流扩散方程的“ 不连续” 4DVar变分融合理想试验

4.1 线性平流扩散方程及变分融合观测算子介绍

4.1.1 4DVar预报模式— 线性平流扩散方程

平流扩散方程是一个抛物线偏微分方程, 其本质上是无散度和忽略压力梯度的Navier-Stocks方程, 定义为:

$\begin{matrix} \{\begin{matrix} \frac{\partial x (s, t)}{\partial t} + a (s, t) \nabla x (s, t) = θ \nabla^{2} x (s, t) \\ x (s, 0) = x_{0} (s) \end{matrix} \end{matrix}$ (12)

式中:a(s, t)表示速度, θ ≥ 0表示黏性常数。简化的线性形式(a为常数)和非黏性(θ =0)方程常被用于资料同化和数据融合新算法的数值模拟中^[29]。本文只考虑其线性形式。基于线性平流扩散和叠加原理, 公式(12)在时间t的解是通过卷积高斯核得到, 即:

$\begin{matrix} D (s, t) = {(4 πθt)}^{- 1 / 2} \exp (\frac{- | s |^{2}}{4 θt}) \end{matrix}$ （13）

式中:标准偏差为 $\begin{matrix} \sqrt[]{2 θt} \end{matrix}$ 。线性幅度变化可以通过使用克罗内克(Kronecker)卷积初始条件得到:

$\begin{matrix} A (s - at) = \{\begin{matrix} 1, s = at \\ 0, otherwise \end{matrix} \end{matrix}$ （14）

文中离散平流扩散模型可以表示为扩散系数(D_0,_i)和初始空间的线性变化幅度(A_0,_i)的积, D_0,_i为托普利兹矩阵^[30]。该试验预报方程为x_i=M_0,_ix₀, 可进一步表示为:

$\begin{matrix} x_{i} = A_{0, i} D_{0, i} x_{0} (15) \end{matrix}$

4.1.2 变分融合观测算子介绍

H= $\begin{matrix} \frac{1}{8} \end{matrix} \begin{matrix} [\begin{matrix} 11111111 & 00000000 & \dots & 00000000 \\ 00000000 & 11111111 & \dots & 00000000 \\ ⋮ & ⋮ & ⋮ \\ 00000000 & 00000000 & \dots & 11111111 \end{matrix}] \end{matrix}$ ∈ R^{n× m}(16)

4.2 变分融合误差协方差和观测值介绍

4.2.1 观测和背景误差协方差介绍

观测和背景误差是资料变分融合重要组成部分, 其决定了系统极小化迭代时的稳定性、解的收敛性和分析场状态质量。对于文中观测误差R, 只考虑固定白噪声的高斯分布, 即R= $\begin{matrix} σ_{r}^{2} \end{matrix}$ I, 其中, v~N(0, R)。

背景误差表示状态变量之间的相关性, 在资料变分融合时, 作为背景状态对目标泛函的权值。一阶自回归(Auto-Regressive, AR(1))高斯马尔可夫过程具有空间相关性的数学模型。本试验将AR作为背景误差; 并在Gaspari等^[31]研究的基础上构建误差协方差模型。定义AR(1)的指数协方差函数为ρ (ϑ )∝ e^{-α |ϑ |}。其中ϑ 表示时间或空间的滞后, 参数α 决定变量间相关性的衰变率, 1/α 通常被称为特征相关长度。该协方差只是函数滞后的大小相关。因离散背景误差协方差是一个埃尔米特矩阵, 可以分解为一个标量标准差和相关矩阵的乘积^[32], 即B= $\begin{matrix} σ_{b}^{2} \end{matrix}$ C_b。

本试验中令σ _b=0.20, σ _r=0.16,

C_b= $\begin{matrix} [\begin{matrix} ρ (0) & ρ (1) & \dots & ρ (m) \\ ρ (1) & ρ (0) & \dots & ⋮ \\ ⋮ & \dots & ρ (1) \\ ρ (m) & \dots & ρ (1) & ρ (0) \end{matrix}] \end{matrix}$ ∈ R^{m× m}。

4.2.2 背景和观测值介绍

为验证本文使用的L1范数正则化的可行性, 共设计了4组对比试验。“ 真实” 状态(即“ 真值” )分别为分段函数(Piecewise Function, PF)、二次函数(Quadratic Function, QF)、正弦函数(Sine Function, SF)和平方指数函数(Square Exponential Function, SEF)。所有初始状态都假定在R²⁰⁴⁸空间, L=2 048; 平流扩散随时间演化的粘度系数θ =4[L/T]; 速度a=1[L/T]。假定变分融合时间间隔介于0至T=500[T]之间, 在此区间内每100[T]步进行观测采样(图1)。

从图1可以看出, 在构建的“ 真实状态” 中有较多点的左右导数不相等, 即在某点x₀有f'( $\begin{matrix} x_{0}^{+} \end{matrix}$ )≠ f'( $\begin{matrix} x_{0}^{-} \end{matrix}$ )。

因变分融合分析必须包括背景场和实际观测值^[4], 且考虑到资料的“ 不连续性” , 文中构建“ 真实状态” 的4种函数设为 $\begin{matrix} x_{0}^{t} \end{matrix}$ , 并加入随机扰动得到“ 背景状态” $\begin{matrix} x_{0}^{b} \end{matrix}$ :

$\begin{matrix} x_{0}^{b} = x_{0}^{t} + \sum_{i} ε_{i} ζ_{i}^{\frac{1}{2}} L_{i} \end{matrix}$ （17）

式中: $\begin{matrix} x_{0}^{b} \end{matrix}$ 为扰动后的背景状态; $\begin{matrix} x_{0}^{t} \end{matrix}$ 为真实状态; L_i和ζ _i分别为背景误差协方差矩阵B的特征向量和特征值; ε _i是均值为0方差为1的高斯随机数。

	Figure Option View Download New Window
	图1 构建“ 真实状态” (a)分段函数(记为PF); (b)二次函数(记为QF); (c)正弦函数(记为SF); (d)平方指数函数(记为SEF)Fig.1 Building true state (a) Piecewise Function (PF); (b) Quadratic Function (QF); (c) Sine Function (SF); (d) Square Exponential Function (SEF)

采用观测算子H(· )和平流扩散方程M(· )对“ 真实状态” 计算得到“ 真实观测值” , 并在此基础上加入随机扰动后得到的值作为实际观测值, 即:

y=H(M( $\begin{matrix} x_{0}^{t} \end{matrix}$ ))+ε _oR^1/2 (18)

式中:y为扰动后的观测值, 即为实际观测值 $\begin{matrix} y_{0}^{o} \end{matrix}$ ; H(M( $\begin{matrix} x_{0}^{t} \end{matrix}$ ))为真实观测值; ε _o为高斯随机误差; R为观测误差。

4.3 “ 不连续” 变分融合理想试验

利用4DVar和R4D-Var 2种方法试验4组独立个例变分融合效果, 结果见图2。其中4DVar和R4D-Var分别为经典变分融合法和经典变分融合耦合L1范数正则项得到的融合结果。

从图2可以看出耦合L1范数正则项(R4D-Var)的融合效果明显优于经典变分(4DVar)融合方法, 尤其对于一些不连续点的融合。分析其原因:①在正则项中分析误差被“ 适当” 地抑制和过滤, 初始状态特征被“ 较好” 地保留; ②经典4DVar一般“ 吻合” 和遵循背景状态而不是提取真实状态; ③经典4DVar在处理白噪声的误差协方差时具有较弱的过滤效果, 存在的冗余信息导致过度拟合; ④经典4DVar基于二次泛函易造成一些“ 拐点” 平滑(拐点即为左右导数不相等的点)^{[4, 13]}。

需要说明的是, 该试验中采用新的融合模型在“ 小波空间” 完成资料融合后, 再采用小波逆变换映射回“ 观测空间” ^[33]。

为验证R4D-Var的融合效果, 本文进一步对100个独立个例开展4组试验, 并采用相对均方误差(Mean Squared Error, MSE_r)、相对绝对误差均值(Relative Mean Absolute Error, MAE_r)和相对偏差(BIAS_r)3种方法进行度量, 结果见表1。3种度量方法分别定义如下^[13]:

MSE_r=‖ $\begin{matrix} x_{0}^{t} \end{matrix}$ - $\begin{matrix} x_{0}^{a} \end{matrix}$ ‖ ₂/‖ $\begin{matrix} x_{0}^{t} \end{matrix}$ ‖ ₂

MAE_r=‖ $\begin{matrix} x_{0}^{t} \end{matrix}$ - $\begin{matrix} x_{0}^{a} \end{matrix}$ ‖ ₁/‖ $\begin{matrix} x_{0}^{t} \end{matrix}$ ‖ ₁

BIAS_r=‖ $\begin{matrix} \overset{{(x}_{0}^{t}}{} \end{matrix}$ ‖ /‖ $\begin{matrix} {\bar{x}}_{0}^{t} \end{matrix}$ ‖ (19)

式中: $\begin{matrix} x_{0}^{t} \end{matrix}$ 表示真实状态, $\begin{matrix} x_{0}^{a} \end{matrix}$ 表示变分融合后的分析状态, $\begin{matrix} \bar{(\cdot)} \end{matrix}$ 表示均值。

	Figure Option View Download New Window
	图2 经典变分融合和耦合L1范数正则化的融合效果对比分析Fig.2 Fusion effect analysis of classic variational fusion and coupling L1 norm regularization method

表1 经典变分和耦合L1范数正则化融合的MSE_r, MAE_r和BIAS_r对比分析 Table 1 Comparative analysis fusion result of classical variational and coupling L1 norm regularization based on the method of MSE_r, MAE_r and BIAS_r

从表1基于一些传统“ 保真度” 度量指标可以看出本文使用的R4D-Var方法整体优于经典4DVar。综合图2和表1结果可以看出, R4D-Var可行, 为不连续资料融合奠定了一定的理论基础。

5 基于“ 不连续” 降水资料3DVar变分融合试验

由于降水资料具有“ 不连续” 和“ 非高斯” 性, 本文只开展降水资料“ 不连续” 的变分融合研究。

5.1 CMORPH反演降水资料订正介绍

本文“ 不连续” 降水试验数据时间段为2011年5月1日至9月30日的CMORPH和安徽省81个国家站观测的小时降水数据。CMORPH反演降水资料从网站直接下载(网站:ftp:∥ftp.cpc.ncep.noaa.gov/precip /global _CMORPH /30min_8 km), 并累加成逐时降水后使用8 km经纬度网格降水场。因该降水资料具有明显系统偏差, 在资料融合前通过GAMMA函数^[34]拟合的概率密度匹配PDF方法进行系统偏差订正。

图3a为砀山站CMORPH和雨量计匹配图; 横坐标表示降水量阈值(单位:mm/h)。图3b为81个站累积概率(Cumulative Probability )在0.99时的偏差(单位:mm)。偏差定义为雨量计减去CMORPH值。以图3a砀山站为例。给出订正过程:第一步, 给定累积概率值如0.99(图中灰色线), 找到红线CMORPH相应累积概率值0.99对应的降水量2.65 mm/h; 第二步, 由0.99找到蓝线对应的CMORPH值降水量值4.85 mm。此时订正量为2.2 mm。由图3b可知, 安徽南部地区CMORPH值偏小, 尤其黄山(30.144° N, 118.16° E)地区。安徽北部CMORPH值偏大, 安徽中部相当。分析其原因与CMORPH得到面降水量有关, 也可能与安徽每年发生的台风雨和梅雨等降水有关。从气候诊断出发, 安徽南部偏湿易发生洪涝, 安徽北部偏干易发生干旱, 安徽中部处的方差变化较小。说明CMORPH资料具有较好的适用性。

5.2 多源降水资料融合试验

进一步对比分析经典3DVar和R3D-Var融合方法的效果(图4)。此时R3D-Var为当M=I时的R4D-Var。图4给出了区域融合站点分布(不同颜色代表站点海拔高度, 颜色超暖则海拔越高)、雨量计观测值(单位:mm/h)分布、CMORPH反演降水值(单位:mm/h)和R3D-Var不同正则化参数(λ )得到的融合场(单位:mm/h)结果。需要说明的是, 此处参考场为中国气象局国家气象信息中心融合产品, 其为CMORPH和地面雨量计观测资料二源融合产品^[14]。地面观测降水资料由国家气象信息中心数据室提供。在操作过程中把分辨率0.1° × 0.1° 的参考场插值到CMORPH和融合场的分辨率, 即0.0625° × 0.0625° 。

从图4可以看出, 雨量计观测给出了单个站的降水分布情况, 无法把握降水的空间分布结构; CMORPH降水资料由红外卫星资料等反演得到, 能够把握降水场的空间结构。CMORPH反演降水和地面雨量计观测值的降水值量级相当, 与R3D-Var融合结果有较好的一致性。并与参考场的形态接近。参考场量级偏小, 其原因与所使用的融合方法(最优插值, OI)有关。区别于参考场更“ 信赖” 地面观测, 从图4e和f看出R3D-Var则较多依赖于背景场CMORPH资料, 且随着正则参数λ 变大与参考场量级更接近, 并且逐步把背景场信息融入到融合场, 能够把握住背景场的结构分布特征。

本文仅对新算法的应用进行研究, 为揭示背景场信息如何作为正则项引入到融合场, 从降水场结构出发, 采用结构相似性(Structural SIMilarity, SSIM)进行度量^[35]。SSIM为局部尺度度量法。对于需要比较的2个数据集U和V的SSIM定义为:

SSIM(u, v)= $\begin{matrix} \frac{(2 \bar{u} \bar{v} + c_{1}) (2 \overset{u}{(} + c_{2})}{({\bar{u}}^{2} + {\bar{v}}^{2} + c_{1}) (\bar{{(u - \overset{̅}{u})}^{2}} + \bar{{(v - \overset{̅}{v})}^{2}} + c_{2})} \end{matrix}$ (20)

式中:c₁和c₂为常数, 用于稳定计算; u和v分别为数据集U和V的部分数据集; $\begin{matrix} \bar{u} \end{matrix}$ 和 $\begin{matrix} \bar{v} \end{matrix}$ 分别为u和v某点周围所选的邻域(如5× 5和7× 7等)均值^{[13, 35]}。

选取28° ~34° N, 104° ~107° E和32° ~38° N, 116° ~119° E 2个区域分析正则化参数λ 值的变化对结果的影响(图5), 即图4中的红色框区域。图5分别给出了λ 值为0.01, 0.02, 0.04, 0.06, 0.08, 0.1, 0.15和0.18的融合结果。图5中标记“ 融合场和CMORPH” 和“ 融合场和参考场” 为R3D-Var得到的融合场分别与CMORPH和参考场的结构相似性, 横坐标中“ 参考场” 表示参考场与CMORPH的结构相似性SSIM(虚线), 该值为比较的“ 准则” 。

从图5中可以看出, 区域1参考场与CMORPH反演降水资料的结构相似性SSIM为0.540。随着正则化参数λ 值增大, 结果与CMORPH和参考场的SSIM值也越来越大, 超过了0.8, 说明文中使用的方法与参考场具有较好的相似性。区域2参考场与CMORPH反演降水资料的结构相似性SSIM为0.436。随着正则化参数变大, 融合场与CMORPH的SSIM越来越大, 说明了正则项能够通过λ 的变化逐渐把背景场信息引入到融合模型中。随着λ 变大, 与参考场的SSIM逐渐变小。结合图4分析可知, 区域2融合降水场具有较好的连续性, 有一条明显“ 雨带” , 而参考场雨带呈现了断裂, 说明文中的方法能够把降水场的结构信息融入到模型中。

	Figure Option View Download New Window
	图3 代表站累积概率匹配图和安徽国家站偏差分布图 (a)CMORPH和雨量计累积概率匹配; (b)CMORPH降水资料对应的国家站偏差值Fig.3 Cumulative probability matching of representative station and bias distribution of national station in Anhui Province (a) Cumulative probability matching of CMORPH and rain gauge; (b) Bias value of CMORPH of corresponding station

	Figure Option View Download New Window
	图4 不同正则化参数的融合降水场与参考场的降水分布结构对比Fig.4 Comparison structure of precipitation fusion based on different regularization parameter and reference field

	Figure Option View Download New Window
	图5 参数不同值融合场与参考场和CMORPH的结构相似性SSIM对比分析Fig.5 Comparative analysis of SSIM between fusion field based on different value of parameters and reference field and CMORPH

6 总结与展望

经典3DVar/4DVar变分融合要求目标代价泛函为凸函数且可微光滑。当资料具有“ 不连续性” 时, 代价函数非光滑。此时经典变分融合基于梯度进行极小化迭代的方法不再适用。文中针对不连续资料进行了变分融合算法的应用研究, 主要结论如下:

(1) 文中将背景资料特点(如不连续性)作为L1范数正则项耦合到3DVar/4DVar模型中。通过线性平流扩散方程作为预报模式, 平滑滤波和采样作为观测算子, 构建4个存在“ 拐点” 的试验, 得到文中使用的方法对这些点的融合具有较好的“ 保真性” 。经典4DVar在这些点处出现了“ Gibbs现象” 或者过度平滑了这些点的信息。由此试验得到此方法在融合一些存在天气现象的“ 离群点” 具有较好的效果。

(2) 在理论分析基础上开展了4组独立试验, 并采用相对均方误差、相对绝对误差均值和相对偏差进行度量, 得到文中使用的方法整体优于经典4DVar。

(3) 在多源降水资料融合时, 基于GAMMA拟合函数PDF法进行CMORPH反演降水资料订正。对比分析得到CMORPH反演降水资料整体比台站观测偏小。

(4) 将经典变分融合耦合L1范数正则项方法用于CMORPH反演降水和台站观测降水融合中。通过结构相似性SSIM度量不同正则化参数λ 值融合结果与参考场、背景场CMORPH的结构相似性, 得到该方法能较好地把背景场中的降水分布结构信息引入到模型, 并表现在融合场中。其融合结果与被融合资料数值量级相当, 而参考场量级偏小。

文中正则项的引入能够逐步引入背景场信息, 能够把握住水汽场的结构, 不仅为一些不连续“ 快变” 量融合, 也为一些连续“ 慢变” 量融合, 如气温、海洋表面温度等融合的研究奠定了理论基础。针对降水资料的非高斯和不连续性融合难题, 尤其是0值的处理, 小波域稀疏性正则化变分融合有望能够解此难点^[36]。经典变分融合一般调用有限内存L-BFGS(Limited-Memory BFGS)极小化算法进行求解。而当目标泛函非光滑时, 文中提到的用于解决大规模非光滑凸函数优化问题的内点算法、近端加速梯度方法和梯度投影等算法具有较好的应用前景。因本文仅对算法的应用进行研究, 且开展的个例研究较少, 代表性不够。另正则化参数也是结合经验给定。后期将开展自适应正则化参数和文中算法大量实际个例研究, 将成果业务化到降水监测、预报和预警中, 以期更好地开展公共气象服务。

致谢:在此衷心感谢Ebtehaj A M博士给予的理论解释和指导, 及评审专家在百忙之中为本文提出的宝贵意见。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	Ghil M, Malanotte-Rizzoli P. Data assimilation in meteorology and oceanography[J]. Advances in Geophysics, 1991, 33(33): 141-266. [本文引用:1]
[2]	Daley R. Atmospheric Data Analysis[M]. Cambridge: Cambridge University Press, 1993. [本文引用:1]
[3]	Kalnay E. Atmospheric Modeling, Data Assimilation and Predictability[M]. Cambridge: Cambridge University Press, 2003: 341. [本文引用:1]
[4]	Wang G, Zhang J W. Generalised variational assimilation of cloud-affected brightness temperature using simulated hyper-spectral atmospheric infrared sounder data[J]. Advances in Space Research, 2014, 54(1): 49-58. [本文引用:6]
[5]	Evensen G. Data Assimilation: The Ensemble Kalman Filter[M]. Berlin/Heiderleng: Springer Science & Business Media, 2009: 307. [本文引用:1]
[6]	Han Pei, Shu Hong, Xu Jianhui. A comparative study of background error covariance localization in EnKF data assimilation[J]. Advances in Earth Science, 2014, 29(10): 1 175-1 185. [韩培, 舒红, 许剑辉. EnKF同化的背景误差协方差矩阵局地化对比研究[J]. 地球科学进展, 2014, 29(10): 1 175-1 185. ] [本文引用:1]
[7]	Talagrand O, Courtier P. Variational assimilation of meteorological observations with the adjoint vorticity equation. I: Theory[J]. Quarterly Journal of the Royal Meteorological Society, 1987, 113(478): 1 311-1 328. [本文引用:1]
[8]	Lorenc A C, Ballard S P, Bell R S, et al. The Met. Office global three-dimensional variational data assimilation scheme[J]. Quarterly Journal of the Royal Meteorological Society, 2000, 126(570): 2 991-3 012. [本文引用:1]
[9]	Law K J H, Stuart A M. Evaluating data assimilation algorithms[J]. Monthly Weather Revew, 2012, 140(11): 3 757-3 782. [本文引用:1]
[10]	Zhang Hua, Xue Jishan, Zhuang Shiyu, et al. Idea experiments of GRAPES three-dimensional variational data assimilation system[J]. Acta Meteorological Sinica, 2004, 62(1): 31-41. [张华, 薛纪善, 庄世宇, 等. GRAPES三维变分同化系统的理想试验[J]. 气象学报, 2004, 62(1): 31-41. ] [本文引用:1]
[11]	Ebtehaj A M, Foufoula-Georgiou E, Zhang S Q, et al. Non-Smooth Variational Data Assimilation with Sparse Priors[J/OL]. 2012. https: ∥www. researchgate. net/ publication/228096140. [本文引用:1]
[12]	Freitag M A, Nichols N K, Budd C J. Resolution of sharp fronts in the presence of model error in variational data assimilation[J]. Quarterly Journal of the Royal Meteorological Society, 2013, 139(672): 742-757. [本文引用:1]
[13]	Ebtehaj A M, Zupanski M, Lerman G, et al. Variational data assimilation via sparse regularization[J]. Tellus A, 2014, 66: 21789, http://dx.doi.org/10.3402/tellusa.v66.21789. [本文引用:12]
[14]	Pan Yang, Shen Yan, Yu Jingjing, et al. Analysis of the combined gauge-satellite hourly precipitation over China based on the OI technique[J]. Acta Meteorologica Sinica, 2012, 70(6): 1 381-1 389. [潘旸, 沈艳, 宇婧婧, 等. 基于最优插值方法分析的中国区域地面观测与卫星反演逐时降水融合试验[J]. 气象学报, 2012, 70(6): 1 381-1 389. ] [本文引用:4]
[15]	Shen Y, Zhao P, Pan Y, et al. A high spatiotemporal gauge-satellite merged precipitation analysis over China[J]. Journal of Geophysical Research: Atmospheres, 2014, 119(6): 3 063-3 075. [本文引用:1]
[16]	Xie P, Xiong A Y. A conceptual model for constructing high-resolution gauge-satellite merged precipitation analyses[J]. Journal of Geophysical Research: Atmospheres, 2011, 116(D21), doi: 10.1029/2011JD016118. [本文引用:1]
[17]	Pan Yang, Shen Yan, Yu Jingjing, et al. An experiment of high-resolution gauge-radar-satellite combined precipitation retrieval based on the Bayesian merging method[J]. Acta Meteorologica Sinica, 2015, 73(1): 177-186. [潘旸, 沈艳, 宇婧婧, 等. 基于贝叶斯融合方法的高分辨率地面—卫星—雷达三源降水融合试验[J]. 气象学报, 2015, 73(1): 177-186. ] [本文引用:2]
[18]	Tikhonov A N, Arsenin V Y. Solutions of Ill-posed Problems[M]. Washington DC: Winston & Sons, 1977. [本文引用:1]
[19]	Wang Gen, Tang Fei, Liu Xiaobei, et al. Application of M-estimators method on FY3B/IRAS channel brightness temperature generalized variational assimilation[J]. Journal of Remote Sensing, 2017, 21(1): 52-61. [王根, 唐飞, 刘晓蓓, 等. M-估计法广义变分同化FY3B/IRAS通道亮温[J]. 遥感学报, 2017, 21(1): 52-61. ] [本文引用:3]
[20]	Miao Chunsheng, Cheng Yuan, Wang Jianhong, et al. Data fusion of offshore SST from China FY and HY2 satellites and its application[J]. Advances in Earth Science, 2015, 30(10): 1 127-1 143. [苗春生, 程远, 王坚红, 等. 中国风云卫星与海洋卫星近海SST资料融合技术及应用研究[J]. 地球科学进展, 2015, 30(10): 1 127-1 143. ] [本文引用:1]
[21]	Rawlins F, Ballard S P, Bovis K J, et al. The Met Office global four-dimensional variational data assimilation scheme[J]. Quarterly Journal of the Royal Meteorological Society, 2007, 133(623): 347-362. [本文引用:1]
[22]	Elad M, Milanfar P, Rubinstein R. Analysis versus synthesis in signal priors[J]. Inverse Problems, 2007, 23(3): 947. [本文引用:1]
[23]	Goldfarb D, Yin W. Second-order cone programming methods for total variation-based image restoration[J]. SIAM Journal on Scientific Computing, 2005, 27(2): 622-645. [本文引用:1]
[24]	Kim S J, Koh K, Lustig M, et al. An interior-point method for large-scale-regularized least squares[J]. IEEE Journal of Selected Topics in Signal Processing, 2007, 1(4): 606-617. [本文引用:1]
[25]	Beck A, Teboulle M. Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[J]. IEEE Transactions on Image Processing, 2009, 18(11): 2 419-2 434. [本文引用:1]
[26]	Nowak R D, Wright S J. Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems[J]. IEEE Journal of Selected Topics in Signal Processing, 2007, 1(4): 586-597. [本文引用:1]
[27]	Hansen P C, O’Leary D P. The use of the L-curve in the regularization of discrete ill-posed problems[J]. SIAM Journal on Scientific Computing, 1993, 14(6): 1 487-1 503. [本文引用:2]
[28]	Nadarajah S. A generalized normal distribution[J]. Journal of Applied Statistics, 2005, 32(7): 685-694. [本文引用:1]
[29]	Smith K S, Marshall J. Evidence for enhanced eddy mixing at middepth in the Southern Ocean[J]. Journal of Physical Oceanography, 2009, 39(1): 50-69. [本文引用:1]
[30]	Chan R H F, Jin X Q. An Introduction to Iterative Toeplitz Solvers[M]. Philadelphia: Society for Industrial and Applied Mathematics, 2007. [本文引用:1]
[31]	Gaspari G, Cohn S E. Construction of correlation functions in two and three dimensions[J]. Quarterly Journal of the Royal Meteorological Society, 1999, 125(554): 723-757. [本文引用:1]
[32]	Williams C K I, Rasmussen C E. Gaussian Processes for Machine Learning[M]. London: The MIT Press, 2006. [本文引用:1]
[33]	Mallat S G. A theory for multiresolution signal decomposition: The wavelet representation[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1989, 11(7): 674-693. [本文引用:1]
[34]	Wang Lei, Chen Rensheng, Song Yaoxuan. Study of statistical characteristics of wet season hourly rainfall at Hulu watershed with Γ function in Qilian Mountains[J]. Advances in Earth Science, 2016, 31(8): 840-848. [王磊, 陈仁升, 宋耀选. 基于Γ 函数的祁连山葫芦沟流域湿季小时降水统计特征[J]. 地球科学进展, 2016, 31(8): 840-848. ] [本文引用:1]
[35]	Wang Z, Bovik A C, Sheikh H R, et al. Image quality assessment: From error visibility to structural similarity[J]. IEEE Transactions on Image Processing, 2004, 13(4): 600-612. [本文引用:2]
[36]	Ebtehaj A M, Foufoula G E, Lerman G. Sparse regularization for precipitation downscaling[J]. Journal of Geophysical Research: Atmospheres, 2012, 117(D8), doi: 10.1029/2011JD017057. [本文引用:1]

1991

0.0

... 气象预测能力、气候和海洋系统建模等方面发挥了核心作用^[1] ...

1993

0.0

... 其基本思想是不断融入状态变量(如风、温度和压力等)观测值到物理状态中(如云和降水等微物理过程),随时间进行迭代通过变分融合系统减少状态变量的不确定性^[2] ...

2003

0.0

... 后验状态(分析场),并用于预测下一个时间的物理状态^[3] ...

2014

0.0

... 虽然各中心在建立时设计方案各有不同,但原理上融合方案都要最小化一个平方型的评价函数,即最小平方估计^[4] ...

... 2 经典3DVar/4DVar变分融合理论分析经典3DVar/4DVar变分融合方法的核心为最小二乘估计理论^[4] ...

... 前期Wang等^[4]把稳健性较强的M#cod#x02014 ...

... ,为保留这些点,本文基于前期Wang等^[4]的研究基础结合Ebtehaj等^[13]提出的方法,将较稳健的L1范数耦合到4DVar作为目标泛函的正则项 ...

... 因变分融合分析必须包括背景场和实际观测值^[4],且考虑到资料的#cod#x0201c ...

... 平滑(拐点即为左右导数不相等的点)^[4,13] ...

2009

0.0

... ,常用方法分为两大类:一是群递归滤波方法,其本质是把随时间演化得到的统计特性融合到系统以追踪最优状态(如基于卡尔曼滤波驱动的数据融合方法)^[5~7] ...

2014

0.0

. 2014, 29(10):1 175-1 185

A comparative study of background error covariance localization in EnKF data assimilation

EnKF同化的背景误差协方差矩阵局地化对比研究

Han Pei , Shu Hong , Xu Jianhui.

韩培, 舒红, 许剑辉

... 区别于OI和卡尔曼滤波等方法^[6],变分具有灵活性,主要体现在2个方面:①可以同时加入不同来源的观测资料 ...

1987

0.0

2000

0.0

... 二是变分法,如三维/四维变分(3DVar/4DVar),其本质是依靠模式直接优化每一瞬间所有可用的观测资料^[8,9] ...

2012

0.0

... 二是变分法,如三维/四维变分(3DVar/4DVar),其本质是依靠模式直接优化每一瞬间所有可用的观测资料^[8,9] ...

2004

0.0

. 2004, 62(1):31-41 DOI:doi:10.11676/qxxb2004.004

Idea experiments of GRAPES three-dimensional variational data assimilation system

GRAPES三维变分同化系统的理想试验

Zhang Hua , Xue Jishan , Zhuang Shiyu

张华, 薛纪善, 庄世宇

This paper introduces a new scheme of three-dimensional variational assimilation(GRAPeS-3DVAR), which is suitable for a grid-point model.The domain includes China area on horizontal dimension and 15 mandatory levels on vertical dimension.The scheme adopts stream function,unbalanced potential velocity,unbalanced geopotential heig htand humidity as analysis variables which are independent of each other.The cost function takes the way ofincremental method.To solve the problems that B -1 are il-conditioned and too large to be computed,the background termis preconditioned using a sequence of variable transform,which improves the convergence in the minimization process and avoids to calculate B -1 directly.EOF decomposition is used as vertical component of analysis variable transform to project onto vertical modes and separate each three-dimensional field into two-dimensional field.Recursive filters is used as horizontal component of analysis variable transform. LBF GS method is used on solution of optimization.Simple geostrophic balanced relationship is used.The system has passed verifications strictly and proved to be credibility.In order to test the functions of the system,A series ofidea experiments were made.The results showa reasonable relationship among the multi-variables and a fast convergence in the minimized process.

文中介绍了一种新的、适合格点模式的三维变分同化方案GRAPeS-3DVAR。该方案采用相互独立的流函数、非平衡速度势函数、非平衡位势和水汽作为分析变量。通过变量变换对目标函数进行预调节,不仅避免了直接计算背景误差协方差逆矩阵的困难,而且改善了Hessian矩阵的性状,提高了收敛速度。采用EOF分解方法,将三维分析变量投影到垂直摸态上,分解成为二维场;水平方向采用数字(递归)滤波器代替矩阵运算,实现和简化了方案的求解。此外,还考虑了质量场和风场之间的平衡约束关系。理想试验结果表明,GRAPeS-3DVAR能够正确地反映多变量之间相互作用关系,收敛迅速,分析结果合理。

... 优化问题,其解决方案是通过定义一个最小加权欧氏距离,即观测和背景与真值的差值估计^[10],权重由背景和观测误差协方差矩阵给定 ...

2012

0.0

... 从估计理论出发,该问题等价于高斯噪声状态下的最大似然(Maximum Likelihood,ML)估计^[11] ...

2013

0.0

... ^[12] ...

2014

0.0

... 资料融合?为此,本文在Ebtehaj等^[13]研究的基础上,显式地把分析状态底层结构信息,即初始背景场作为正则项耦合到经典变分融合框架中 ...

... 本文基于Ebtehaj等^[13],Tikhonov等^[18]和王根等^[19]的研究方法得到新的变分融合法,为#cod#x0201c ...

... 资料底层分析状态的先验知识作为L1范数正则项^[13],通过线性平流扩散方程进行理想试验(同化的背景资料和观测资料不连续,且具有强的#cod#x0201c ...

... 则观测与状态的关系为公式(1)^[13]: ...

... 3DVar/4DVar变分融合归结为在初始时刻极小化以下目标泛函^[13,22]: ...

... 根据最大后验(Maximum A Posteriori,MAP)估计和贝叶斯定理,在4DVar代价函数中加入正则化约束项,则有^[13,23]: ...

... 本文中基于梯度投影算法(Gradient Projection,GP)进行迭代求解^[13,26] ...

... 基于拉格朗日乘子法公式(8)的约束问题转化为求解无约束问题,则有^[13]: ...

... ,为保留这些点,本文基于前期Wang等^[4]的研究基础结合Ebtehaj等^[13]提出的方法,将较稳健的L1范数耦合到4DVar作为目标泛函的正则项 ...

... 平滑(拐点即为左右导数不相等的点)^[4,13] ...

... 3种度量方法分别定义如下^[13]: ...

... 7等)均值^[13,35] ...

2012

0.0

. 2012, 70(6):1 381-1 389 DOI:doi:10.11676/qxxb2012.116

Analysis of the combined gauge-satellite hourly precipitation over China based on the OI technique

基于最优插值方法分析的中国区域地面观测与卫星反演逐时降水融合试验

Pan Yang , Shen Yan , Yu Jingjing

潘旸, 沈艳, 宇婧婧

In order to develop a high-quality Chinese precipitation product at the highly spatial temporal resolution, the optimum interpolation (OI) technique was adopted to combine the CMORPH with the Chinese Precipitation Analyses (CPA)，which was based on the hourly gauged rainfall from about 30,000 automatic weather stations (AWS), at 0.1°lat/lon resolution. The CMORPH was set to the first guess, while the CPA was used as observations to modify the first guess. The error statistics of the first guess and observations were defined based on the data from June to August of 2009 and the error structures were employed to compute the weights. The examination of the combined precipitation (COMB) at 00:00 UTC 2 July 2009 showed that the CMORPH at the target grid was able to be improved by the CPA if the gauge within a searching radius was available. Otherwise, the CMORPH was remained where there was no gauge within the searching radius. The independent validation result of the COMB from June to August of 2009 was that, the Bias/RMSE/RE of the COMB were -0.004 mm/h, 1.271 mm/h and 15.964%, respectively, with the averaged spacial correlation coefficient of 0.778. It was indicated that the COMB had a better quality than the CMORPH, as well as the FY-2C merged precipitation product.

为了发展一套适用于中国区域的高分辨率（0.1°×0.1°）逐时降水产品，以CMORPH卫星反演降水为背景场，以基于3万个自动气象站观测的逐时降水量分析的中国降水格点分析产品（Chinese Precipitation Analyses，CPA）作为地面观测场，采用最优插值方法对二者进行了融合试验。用2009年6—8月的样本统计分析了卫星反演与地面观测降水的误差及其协相关形式，按照误差结构来分配权重。融合试验的个例检验表明，该方案在有站点的地区能较好地引入地面观测信息，在没有站点观测的地区则保留CMORPH的原始信息，最终形成一套覆盖中国区域的高时空分辨率的降水场。2009年6—8月独立样本检验的统计结果也表明，该融合产品的平均偏差、均方根误差、相对误差分别为-0.004 mm/h、1.271 mm/h和15.964%，平均空间相关系数达到0.778，与融合前CMORPH的各统计值相比，改进幅度基本都超过了50%，且与风云系列卫星的同类型产品相比精度也有一定程度的提高。

... ^[14~16]和#cod#x0201c ...

... 融合研究较多的是融合CMORPH反演降水资料与地面站观测资料,采用概率密度匹配法(Probability Density Function matching method,PDF)和最优插值(Optimal Interpolation, OI)^[14] ...

... )可以假设为线性算子^[14] ...

... 需要说明的是,此处参考场为中国气象局国家气象信息中心融合产品,其为CMORPH和地面雨量计观测资料二源融合产品^[14] ...

2014

0.0

2011

0.0

... ^[14~16]和#cod#x0201c ...

2015

0.0

. 2015, 73(1):177-186 DOI:10.11676/qxxb2015.010

An experiment of high-resolution gauge-radar-satellite combined precipitation retrieval based on the Bayesian merging method

基于贝叶斯融合方法的高分辨率地面#cod#x02014;卫星#cod#x02014;雷达三源降水融合试验

Pan Yang , Shen Yan , Yu Jingjing

潘旸, 沈艳, 宇婧婧

In order to develop a method for combining different precipitation sources such as gauge, radar and satellite with high resolution, the Bayesian merging method, which has been used to retrieve high-resolution surface moisture from radar and satellite sensors, was adopted to combine gauge, radar and satellite precipitation data at hourly and 0.05°×0.05° lat/lon resolutions over the Jiang-huai region. The errors statistics of the three sources of precipitation was performed via the samples in August of 2009. The CMORPH was set to the first guess, and the radar and gauge precipitations were merged as the observations together. The verification at the independent stations showed that the Bayesian merging method can be effectively applied to the combination of gauge-radar-satellite precipitations, and the accuracy of the combined precipitation was higher than any of the three sources of precipitation.

为了探讨一种适用于区域性的地面、雷达、卫星等多源降水资料融合的方法,一种曾用于高分辨雷达、卫星土壤湿度产品反演的贝叶斯融合(Bayesian Merging)方法被尝试应用于江淮地区1 h-0.05°×0.05°经纬度高分辨率的雷达估测降水、卫星反演降水与地面站点观测降水3种资料的融合。在应用该方法时,通过2009年8月样本统计分别估计卫星和雷达反演降水的误差关系,通过曲线拟合建立误差方程,并以卫星资料作为背景场,但在融合时将雷达估测降水作为新的观测信息与地面观测降水同时引入。融合试验检验结果表明:贝叶斯融合方法能够有效实现雷达、地面、卫星3种不同来源资料的融合,该方法生成的多源融合产品的精度均优于任何单一来源的降水产品。

... 融合^[17] ...

... 其具体思路为先采用PDF方法订正雷达和卫星估测降水的系统误差,再采用贝叶斯融合方法将三者进行融合^[17] ...

1977

0.0

... 本文基于Ebtehaj等^[13],Tikhonov等^[18]和王根等^[19]的研究方法得到新的变分融合法,为#cod#x0201c ...

2017

0.0

... 本文基于Ebtehaj等^[13],Tikhonov等^[18]和王根等^[19]的研究方法得到新的变分融合法,为#cod#x0201c ...

... 而对于卫星辐射资料,该算子为辐射传输模式^[19,21] ...

... 通过求解演变模型线性化后的切线性模型的伴随模式,即可得到相应的雅克比矩阵^[19] ...

2015

0.0

. 2015, 30(10):1 127-1 143

Data fusion of offshore SST from China FY and HY2 satellites and its application

中国风云卫星与海洋卫星近海SST资料融合技术及应用研究

Miao Chunsheng , Cheng Yuan , Wang Jianhong

苗春生, 程远, 王坚红

... ,如气温、海洋表面温度( Sea Surface Temperature,SST)^[20]等一些#cod#x0201c ...

2007

0.0

... 而对于卫星辐射资料,该算子为辐射传输模式^[19,21] ...

2007

0.0

... 3DVar/4DVar变分融合归结为在初始时刻极小化以下目标泛函^[13,22]: ...

2005

0.0

... 根据最大后验(Maximum A Posteriori,MAP)估计和贝叶斯定理,在4DVar代价函数中加入正则化约束项,则有^[13,23]: ...

2007

0.0

... 目前国际上用于解决大规模非光滑凸函数的优化问题主要有内点算法(interior point algorithms)^[24]和近端加速梯度方法(accelerated proximal gradient methods)^[25]等 ...

2009

0.0

2007

0.0

... 本文中基于梯度投影算法(Gradient Projection,GP)进行迭代求解^[13,26] ...

1993

0.0

... 非负^[27] ...

... 参数的确定是通过统计交叉验证获得^[27] ...

2005

0.0

... 正则化约束是指状态量关于先验知识的概率分布^[28],即p(x)#cod#x0221d ...

2009

0.0

... =0)方程常被用于资料同化和数据融合新算法的数值模拟中^[29] ...

2007

0.0

... 文中离散平流扩散模型可以表示为扩散系数(D_0,_i)和初始空间的线性变化幅度(A_0,_i)的积,D_0,_i为托普利兹矩阵^[30] ...

1999

0.0

... 并在Gaspari等^[31]研究的基础上构建误差协方差模型 ...

2006

0.0

... 因离散背景误差协方差是一个埃尔米特矩阵,可以分解为一个标量标准差和相关矩阵的乘积^[32],即B=σb2C_b ...

1989

0.0

... ^[33] ...

2016

0.0

. 2016, 31(8):840-848 DOI:10.11867/j.issn.1001-8166.2016.08.0840.

Study of statistical characteristics of wet season hourly rainfall at Hulu watershed with #cod#x00393; function in Qilian Mountains

基于#cod#x00393; 函数的祁连山葫芦沟流域湿季小时降水统计特征

Wang Lei , Chen Rensheng , Song Yaoxuan.

王磊, 陈仁升, 宋耀选

1.Qilian Alpine Ecology and Hydrology Research Station,Cold and Arid Regions Environmental and Engineering Research Institute, CAS, Lanzhou 730000,China; 2.Key Laboratory of Inland Ecohydrology, CAS,LanZhou 730000,China

The probability distribution of the wet season hourly precipitation is the important basis for the study of the precipitation distribution, especially in mountainous areas. Hulu watershed is the study area located in the upper reaches of Heihe River, Qilian Mountains. By adopting the maximum likelihood estimation, the shape parameter α and scale parameter β of 6 stations were obtained with observed wet season (May to September) half hourly data, and different intensity precipitation probability density distribution, cumulative probability density and probability of precipitation and elevation and precipitation relationship were analyzed. The shape parameter α and scale parameter β is significantly negatively correlated, shape parameter α and average hourly precipitation distribution is consistent. Local topography is also an important factor to affect the precipitation redistribution and the probability distribution of precipitation in Hulu watershed. In addition to the increase of precipitation events, the probability of 1~3 mm mm/h precipitation increases with the altitude, which is the main reason for the increase of precipitation with altitude.

高山区降水主要集中在湿季(5~9月),湿季小时降水的概率分布统计特征是研究山区降水分布的重要基础。选择祁连山中段典型流域为研究区,使用葫芦沟流域6个观测点2015年5~9月半小时降水数据,采用极大似然法对影响Γ分布函数的形状参数 α 和尺度参数 β 进行估计,并且对不同强度降水概率密度分布、累计概率密度以及降水概率与海拔和降水量关系进行分析。结果表明:形状参数 α 和尺度参数 β 呈明显负相关,形状参数α与小时平均降水量分布大致相同;在葫芦沟流域,除海拔因素外,局地地形也是影响降水再分配和降水概率分布的重要因素;在该流域,除降水事件增加外, 1~3 mm/h强度降水概率随着海拔增加而增加也是降水量随海拔增加而增加的主要原因。

... 因该降水资料具有明显系统偏差,在资料融合前通过GAMMA函数^[34]拟合的概率密度匹配PDF方法进行系统偏差订正 ...

2004

0.0

... 本文仅对新算法的应用进行研究,为揭示背景场信息如何作为正则项引入到融合场,从降水场结构出发,采用结构相似性(Structural SIMilarity,SSIM)进行度量^[35] ...

... 7等)均值^[13,35] ...

2012

0.0

... 针对降水资料的非高斯和不连续性融合难题,尤其是0值的处理,小波域稀疏性正则化变分融合有望能够解此难点^[36] ...