大气海洋数据同化方法变分法课件

大气/海洋数据同化方法（3）变分法大气/海洋数据同化方法（3）1数据同化的目的资料同化:数值天气预报问题是一个初/边值问题,初始场的精确性直接决定着预报的精确性.通过模式预报和观测的统计结合得到大气初始场的过程称为数据同化.Purposeofdataassimilation:usingalltheavailableinformationtodetermineasaccuratelyaspossiblethestateoftheatmosphericflow.---Talagrand(1997)数据同化的目的资料同化:数值天气预报问题是一个初/边值问题,2两点要注意的问题：（1）何谓“精确的大气状态”？要考虑大气运动的多尺度特征，这里的精确（或真实）依赖于研究对象的特征尺度。（2）所有有效信息包括什麽？观测资料（包括它的误差）和控制大气运动的物理规律两个方面。两点要注意的问题：3发展历史简要回顾早期的客观分析和数据同化方法:Panofsky(1949)：提出的2维全局多项式插值方法是第一个客观分析方法；GilchristandCressman(1954)：提出局地多项式插值方法；BergthorssonandDoos(1955)：指出客观分析中应该给出所有格点的初猜值来弥补观测的不足，由此发展了逐步订正法（Cressman1959,Barnes1964,1978）。逐步订正法采用了短期预报的结果作为初猜值，又不断插入观测（6小时一次），这样的循环过程就构成四维资料同化；Gandin(1963)：提出最优插值方法；Sasaki(1970)：首先将变分法引入初始化过程；HokeandAnthes(1976)：提出Nudging（牛顿张弛）四维变分方法。发展历史简要回顾早期的客观分析和数据同化方法:4LewisandDerber(1985),LeDimetandTalagrand(1986),CourtierandTalagrand(1987)提出了四维变分同化方法(4DVAR),随后ECMWF和NMC相继在业务上采用三维和四维变分同化方法。四维变分同化方法是一种最优控制方法，三维变分同化（3DVAR)和最优控制方法等价，是一种最优估计方法。LewisandDerber(1985),LeDimet5大气海洋数据同化方法变分法课件6

变分方法在大气与海洋领域的发展历史20世纪70年代到80年代，在Marchuk及其同事(1975)提出在气象学中运用伴随方程的思想和在Sasaki(1958,1969,1970)提出在气象学中运用变分方法的思想基础之上，Penenko(1976)，LeDmiet(1986)和Talagrand(1987)分别提出了用动力模式作为约束条件构造变分问题(强变分约束)，并用伴随方程去求其叠代解的数据同化新思路，这就是四维变分（4DVAR）伴随同化技术（变分伴随法）。变分伴随法是变分原理和最优控制论（方差最小化）相结合产生的一种方法。

变分方法在大气与海洋领域的发展历史7变分问题的起源问题是数学发展的源泉！著名的“最速降线”问题（TheBrachistochroneProblem）约翰·伯努利（JohannBernoulli，1667－1748）1696年向全欧洲数学家挑战，提出一个难题：“设在垂直平面内有任意两点，一个质点受地心引力的作用，自较高点下滑至较低点，不计摩擦，问沿着什么曲线下滑，时间最短？”变分问题的起源问题是数学发展的源泉！著名的“最速降线”问题（8给出答案的人1:罗比塔（GuillaumeFrancoisAntoniedel‘Hospital,1661-1704）2:雅可比·伯努利（JacobBernoulli,1654-1705）3:莱布尼茨（GottfriedWilhelmLeibniz,1646-1716）4:牛顿（IsaacNewton,1642—1727）约翰的解法比较漂亮。雅可布的解法虽然麻烦与费劲，却更为一般化。后来欧拉（EulerLonhard，1707～1783）和拉格朗日(Lagrange,JosephLouis，1736-1813)发明了这一类问题的普遍解法，从而确立了数学的一个新分支——变分学。给出答案的人1:罗比塔（GuillaumeFrancoi9

瑞士数学家欧拉1744年出版《寻求极大或极小性质的曲线的技巧》，由此

变分法作为一个新的数学分支诞生。这本书给欧拉带来声誉，一度他被视为当时“活着的最伟大的数学家”。瑞士数学家欧拉1744年出版《寻求极大或极小性质的曲线10问题之23：

变分法的进一步开展。

1900年，德国大数学家

D.

Hilbert在巴黎举行的国际数学会议上，发表了〈数学问题〉的专题演讲，提出了23个问题。

问题之23：

变分法的进一步开展。11泛函介绍泛函介绍12大气海洋数据同化方法变分法课件13大气海洋数据同化方法变分法课件14大气海洋数据同化方法变分法课件15请用变分法解

“最速降线”问题！请用变分法解

“最速降线”问题！16满足条件：满足条件：17约束问题：

逐步订正松弛逼近Kalman滤波变分问题的约束条件是模式，或是物理关系。约束问题：18问题的提出：让满足，使最小！

问题的提出：让满足，使19怎么解决？途径：通过构造Lagrange函数化约束问题为无约束问题。怎么解决？20大气海洋数据同化方法变分法课件21弱变分约束问题的提出（变分调整）：约束条件不是整个模式，而是较简单的物理诊断关系。如地转关系不能再运用伴随方程，因为无时间倾向项！弱变分约束问题的提出（变分调整）：22解决途径：1）构造泛函（目标函数）J(u,v,p)=(u-u’)(u-u’)+(v-v’)(v-v’)+(p-p’)(p-p’)2）考虑约束关系，构造Lagrange函数

E(x,y,u,v,p,px,py)=J(u,v,p)+aFxFx+bFyFy数学问题：求使Emin时的(u,v,p)解决途径：233）通过Euler方法，化变分问题为微分问题我们要求u,v,p3）通过Euler方法，化变分问题为微分问题243）将E(x,y,u,v,p,px,py)代入，得到事先给定两个系数：ab它们可调！3）将E(x,y,u,v,p,px,py)代入，得到事先给定25大气数据的三维变分分析方法

大气数据的三维变分分析方法261三维变分分析

如果已知大气的观测y0，背景场xb，那末按照线性估计理论，在统计意义下x的最佳估计（分析场）是

它是下面目标函数的极小点：

H是由计算y的算子(H(x)=y),称为观测算子（可能是简单的内插算子或复杂的模式）。直接由（a）求x是非常困难的，一般用我们前面讲过的下降算法寻找J的极小点。计算梯度的公式是：

其中，是观测算子的切线性算子。1三维变分分析27Courtier等（1994）为了解决四维变分同化方法耗费计算量过大的问题，提出了一种所谓增量方法。可以将(b)写为扰动量（增量）形式：（c）

这里δx=x-xB,d=H(xB)－y,H’是H的切线性算子（Jacobian矩阵）。(对H(x)在xb作泰勒展开,截取前两项)梯度为（d）不过它只给出近似解，（是第n次迭代的值）要将δx+xB

重新赋给xB进行多次迭代计算，才能够给出更精确的解。因为（c）包含了近似：。在H是线性或者接近线性时，这一近似成立，无需迭代计算。Courtier等（1994）为了解决四维变分同化方法耗费计28三维变分同化的困难及解决办法

一般认为不同的观测之间(包括不同变量间)误差的水平方向不相关,在垂直方向可能相关，所以距阵O是对角阵或者接近对角阵，求逆不困难。三维变分同化的困难及解决办法29三维变分同化的困难及解决办法

但是B距阵就不能看为对角阵，而它是一个维数很高的距阵，求逆实际上是不可能的。这是三维变分同化遇到的一大困难。为了解决这一困难，首先是希望分析变量之间（从而他们的误差）相互独立。这可以使背景误差协方差矩阵成为事实上的块对角矩阵）。通常的做法是取模式变量的非平衡部分作为分析变量。它们是互不相关的。这样做还容易控制非平衡模态的增长。

三维变分同化的困难及解决办法30实际大气中，水平风的两个分量是不独立的，温度和高度（气压）是密切相关的，风和高度也有密切相关。实际大气中，31可以考虑作下面的变换原来的大气变量为:现在变换为：

用r代替q是为了避免q随高度变化剧烈带来的的问题。而真正的分析变量是xp的非平衡部分，设它为：带下标的变量表示变量的非平衡部分，即为速度势和位势高度场与流函数平衡的部分。（可以近似令但也可保留），而可以由平衡方程导出：可以考虑作下面的变换32背景误差协方差矩阵的表达式可以写为块对角矩阵： 是单变量误差协方差矩阵，它可以对预报误差作统计分析给出。NCEP（前身为NMC）利用同一时刻不同时效（24和12小时）的预报之差作为预报误差的近似：

这里下标i

和l表示水平方向与垂直方向的格点标号.背景误差协方差矩阵的表达式可以写为块对角矩阵： 是单变量误差33下面我们讨论几种避免求B的逆的办法（1）

Lorenc(1988),DerberandRosati(1989)建议作变换目标函数成为：

梯度：

如果令x的初猜值为xB，即v的初猜值为零，那末实际并不用计算。按照下降算法迭代给出最后的v后，由x=xB+Bv得到分析场。不过在实际问题中B是一个维数很高的距阵（可达107×107），而且每次迭代都要作Bv的运算，占用空间很大，计算量也不小。下面我们讨论几种避免求B的逆的办法目标函数成为：梯度：如34

(2)运用过滤算子的变分分析方法：

分析矩阵B和向量v相乘的作用。相乘后v的第j个分量是

在均匀格点下，上式相当于一个空间滤波运算，是滤波系数，过滤范围是整个模式区间。按照这样的分析，可以用空间滤波运算代替B，这样的方法简称VAF（variationalanalysisusingafilter）。计算过程如下：初猜值v=0,梯度：这里G是一个空间滤波算子，它可以根据B的特性设计。

35例如在二维空间，认为背景场误差各项均匀同性，可用下面的高斯分布表示背景场误差相关函数（Daley,1991）:

是格点i,j之间的距离，是背景场误差标准差。这样矩阵B和v的乘相当于用一个高斯滤波器作用到v上，实际计算时应用的是一个截断的高斯滤波器，可以大大减少内存和计算量。例如在二维空间，认为背景场误差各项均匀同性，可用下面的高斯分36(3)Lorenc(1997)研究了用递归滤波模拟的可能性。

让作变换，则目标函数成为梯度是：

现在的问题是U的具体形式如何给出。前面我们指出，矩阵与某个向量的乘法可以用滤波运算逼近。如果认为背景场误差相关函数可用形如(8.2.7)的高斯分布表示，考虑一维问题，令F＝Asin(kx)代入高斯滤波器（8.2.5）可以计算对波数为k的波的谱响应函数(振幅响应因子)：(3)Lorenc(1997)研究了用递归滤波模拟的可能性37设计一个递归滤波器定义如下：考虑一格距为的一维格点场，是格点i的原始值，是经过从到I滤波后的值，C是经过左右两个方向滤波后的值，是滤波系数。则递归滤波器定义为：将两式合并，有：(8.2.12)(8.2.12)亦称反滤波器，从该式中可以求得其谱响应函数，并进而得到递归滤波器的谱响应函数（为反滤波器谱响应函数的倒数）。

设计一个递归滤波器定义如下：38N次递归滤波的振幅响应因子为：如，则有：（8.2.16）前面给出的高斯谱响应函数（8.2.11）为：（8.2.11）由于故递归滤波器（8.2.12）的结果需要乘一个因子对（8.2.11）关于作Taylor展开：(8.2.17）（8.2.17）与（8.2.16）对比，可得：