四维变分资料同化

切线性模式的检验: 切线性模式的检验:
φ(α) =
Qr (z +αh) Qr (z)
αP h r
=1+ O(α)
Qr
P r
分别是非线性和切线性模式在r时刻 分别是非线性和切线性模式在 时刻 的预报. 的预报.
伴随模式的检验: 伴随模式的检验: 为输入积分切线性模式y=Prz, 以y为输入反向积 以z为输入积分切线性模式 为输入积分切线性模式 为输入反向积 分伴随模式 w=PrTy. 应该有yTy=zTw, 因为
(4)显式四维变分资料同化方法 (邱崇践 张蕾 ) 邵爱梅,中国科学D,2007.5) 将SVD技术用于四维空间的预报集合, 这样得到的 奇异向量就不但可以表征模式变量的空间结构, 也可以包含其时间演变特征. 它在一定的近似程度内保留了传统的4DVAR的基 本优点,但是实现起来会非常简单. 分析增量按照奇异向量展开:(V是4维空间奇异 向量构成距阵) δx = x x ≈ r α v (t,r) =V α
(2)反向4DVar(王斌 赵颖,气象学报 2005.5) ) 提出了映射观测的新概念和反向四维变分资料同化的 新思路,并以此为基础建立了三维变分映射资料同化 (缩写为3DVM ). 由3DVM得到的初值不在同化窗口的 始端,而在窗口的末端. 3DVM所花计算时间只需4DVar 的1/7. (3)4DVAR和EnKF结合 ) 和 结合
(2.3)的切线性方程 的切线性方程
δu δu u +u +δu = 0 t x x
定义伴随方程: 定义伴随方程:
* *
(2.4) )
(δu ) (uδu ) * u ) +δu = w u(x, t) uobs )] (2.5) [ t x x
δu*和(2.5)乘 δu* 相减 在整个区间 将(2.4)乘 ) ) 积分 δu*δu (uδuδu*)
以 δ~为初值反向积分伴随模式到t=0. 最后 xr
~T J o = ∑ P δ~r x r
R
实际过程由R开始积分到下一个观测时间R-1 得到
o o J 相对于 R-1的梯度 JR,R1 ,从 JR,R1 +δ~R1 相对于x x
r =0
再积分伴随方程到R-2, …… 再积分伴随方程到 ,
伴随程序的书写技巧和检验. (大气模式, 伴随程序的书写技巧和检验. 大气模式, 资料同化和可预报性,气象出版社, 资料同化和可预报性,气象出版社,2005) ) 伴随模式的解析形式只能作为理论推导用, 伴随模式的解析形式只能作为理论推导用,实际 问题是离散化的, 问题是离散化的,预报程序中还有些是不能写成 解析公式的.要保证相应的伴随模式严格成立, 解析公式的.要保证相应的伴随模式严格成立, 通常的作法是先根据原模式计算程序写出切线性 模式程序, 模式程序,再直接根据切线性模式程序一一对应 地写出伴随程序. 地写出伴随程序.一个天气预报模式的程序有上 万条语句,首先写出他的切线性模式程序, 万条语句,首先写出他的切线性模式程序,然后 根据切线性模式程序写出伴随程序, 根据切线性模式程序写出伴随程序,工作量是巨 大的. 大的. 按照一定的规则来写. 按照一定的规则来写.
( f (x), g(x)) = ∫ f (x)g(x)dx
( x, y) = xT y
显然矩阵算子A的伴随算子是 显然矩阵算子 的伴随算子是AT, 的伴随算子是
x Ay = y A x
T T T
计算过程: 计算过程: (1)给出u0初猜值积分模式(2.3)得到u(x,t).计算 ) 初猜值积分模式( ) 目标函数 u u +u = 0 t x (2)从u*(x,T)=0 出发积分伴随方程得到u*(x,0) )
zTPrTy=yTPrz
梯度检验: 梯度检验:
J ( x0 +αh) = J ( x0 ) +αhT J ( x0 ) + O(αh)
φ(α) =
J ( x0 +αh) J ( x0 )
αh J ( x0 )
T
=1+ O(α)
h = J ( x0 ) / J ( x0 )
是单位向量
5 关于 关于4DVAR的评论 的评论 好处:利用完整的模式方程约束,同化场动力上 好处:利用完整的模式方程约束, 协调,可以同时同化多时刻资料, 协调,可以同时同化多时刻资料,容易加 入其它约束.已经成功运用. 入其它约束.已经成功运用. 问题:不考虑模式误差; 问题:不考虑模式误差;
bHale Waihona Puke Baidu
∑
k =1
k k
代入目标函数可以直接求梯度
1 T J (α) = (M 1)αT ∧2 α + ∑[H( xB +V α) y]r R1[H( xB +V α) y]r 2 r
再
见
�
目的:给出大气,海洋,陆面…状态的最好估计 状态的最好估计, 目的:给出大气,海洋,陆面 状态的最好估计, 为预报和分析研究提供必要的数据. 为预报和分析研究提供必要的数据.
为什么要用预报模式? 为什么要用预报模式? (1)观测不足 ) (2)观测有误差 ) 后果:变量间的不协调造成预报的振荡) (后果:变量间的不协调造成预报的振荡) 预报模式给我们提供什么? 预报模式给我们提供什么? (1)模式作出的预报为同化提供初猜场(背景场) )模式作出的预报为同化提供初猜场(背景场) (2)模式在不同点的变量之间以及各个变量之间建 ) 立了联系
1 J= ∫ 20
T
b
w u(x, t) uobs )]2 dxdt ∫ [
a
由(2.7)看到: )看到:
u0 J = δu (x,0)
*
(2.8)
伴随算子的定义: 伴随算子的定义:
(f, Lg )=(g, L*f ), ( f,g)内积 , , )
(1) 函数空间内积 (2) 向量空间内积: 向量空间内积:
J (u)= {w [u(x) uobs (x)]2 + w2[2u / x2 ]2}dx ∫ 1
b a
(2.1) ) J 是u的泛函,依赖于u 在(a,b)区间的所有取值 区间的所有取值. 的泛函, 区间的所有取值 (1阶)变分:δJ = J (u +δu) J (u) 对 δu的线性部分 阶 变分:
T
(1)
= (I KH)B
解(1)等价于极小化下面的目标函数 (cost ) function 代价函数 代价函数) 1 1 T 1 J = ( xB xa ) B ( xB xa ) + ( yo Hxa )T Q1( y0 Hxa ) 2 2 (2) 变分方法: 极小化( ) 变分方法: 极小化(2) 4维变分方法: 极小化(3) 维变分方法: 极小化( ) 维变分方法
(δu*) (uδu*) * u obs +δu = w u(x, t) u )] [ t x x
(3)根据目标函数值和梯度找到新的估计u0 ) (4)重复(1)-( )迭代. )重复( )-(3)迭代. )-(
3 4DVAR实际计算过程 回到离散情况 : 实际计算过程
1 1 R T 1 o o J ( x0 ) = ( x0 xB ) B ( x0 xB ) + ∑ (Hxr yr )T Qr1(Hxr yr ) 2 2 r=0 = J B + Jo
先看连续情况.反演初值的一个例子: 先看连续情况.反演初值的一个例子:目标泛函
1 J= ∫ 20
T
b
w u(x, t) uobs )]2 dxdt ∫ [
a
(2.2) ) (2.3) )
u u +u = 0 t x
找到最优的u0 让(2.2)极小 极小
定解条件: 定解条件: u(x,0) = u0 (x), u(a, t) = f (t), u(b, t) = g(t)
∫∫ ∫∫
t
dtdx + ∫∫ u
obs
x
dxdt =
w u(x, t) u )]δudxdt [
(2.6) )
考虑边界 δu = 0 以及令 δu*(x,T) = 0
得到
∫ δu (x,0)δu(x,0)dx = ∫∫
* a
b
w u(x, t) u )]δudxdt [
obs
= δJ
(2.7)
B 距阵不随时间变化(为了减少求逆困难 距阵不随时间变化(
要简化B); 伴随模式程序编写维护工作量大.
6 对改进 对改进4DVAR所作的努力 所作的努力 (1)增量方法,减少计算量(用低分辨率模式) )增量方法,减少计算量(用低分辨率模式)
1 T 1 1 R o J (δx0 ) = δx0 B δx0 + ∑ (Hδxr + dr )T Qr1(Hδxr + dr ) 2 2 r=0 dr = HxB yr
1 1 K o o J ( x0 ) = ( xB x0 )T B1( xB x0 ) + ∑ ( yk Hxk )T Qk 1( yk Hxk ) 2 2 k =0
(3) ) 预报模式: 预报模式:
xk = Mk x0
(4) )
示意图: 示意图:
如何求极小?下降算法 如何求极小?下降算法.
四维变分资料同化
( 4DVAR, four-dimensional variational data assimilation method )
兰州大学大气科学学院 邱崇践
1,有关资料同化的基本知识
资料同化: 资料同化 在积分描写动力系统演变过程的数学模 预报模式)的同时,不断吸收观测资料, 式(预报模式)的同时,不断吸收观测资料,给出 系统状况的一个估计. 系统状况的一个估计.
(需要梯度) 需要梯度)
最速下降法,共轭梯度法,拟牛顿法 最速下降法,共轭梯度法,拟牛顿法…
得到搜寻方向后成为一维寻优. 得到搜寻方向后成为一维寻优. 多项式逼近: 多项式逼近:
2,变分方法和伴随模式 , 变分方法是求泛函极值的有力工具. 变分方法是求泛函极值的有力工具. 是求泛函极值的有力工具 泛函的一个例子: 泛函的一个例子:
(成为求多元函数极小值 . x J B = B1(x0 xB ) 成为求多元函数极小值). 成为求多元函数极小值
0
看一个时间r 看一个时间 ′ δJro = (x0 Jro )T δx0 = Hr [Or 1(Hr xr - yr )]T δxr 切线性模式: 切线性模式: δxr = P1δxr1 = P1P2δxr2 = ...... = r r r ~ P1P2...P δx0 ≡ P δx0 r r 0 r o T 1 T ~ ′ (x0 Jr ) δx0 = Hr [Or (Hr ( xr ) - yr )] Pδx0 r
~T ′ Jro = Pr HrTOr1(Hr ( xr ) - yr )
~ T T 1 ′ J = Pr Hr Or (Hr ( xr ) - yr ) 是什么? 是什么?
o r
~ δxr = P 1 Pr2 ...P δx0 ≡ P δx0 r 0 r ~T ~ Pr δxr = PT PT ...PrT 1δ~r x ′ δ~r = HrT Or1(Hr ( xr ) - yr ) x 0 1
给出的最好估计是什么? 综合观测和背景场 给出的最好估计是什么? 背景场 xB, 观测场yo,分析场 a. 最大似然估计: 观测场 分析场x 最大似然估计: 分析场
xa = xB + BH (HBH + Q) ( yo HxB )
T T
1
= xB + K( yo HxB )
(卡尔曼滤波) (条件?) B: 背景场误差协方差矩阵;Q: 观测场误差协方 背景场误差协方差矩阵; 差矩阵; 观测算子( 差矩阵;H: 观测算子(Hx=y) 误差协方差矩阵: 误差协方差矩阵 B=<ε εT>, bi,j=< εi εj> 分析场的误差协方差矩阵 误差协方差矩阵: 分析场的误差协方差矩阵 P = εaεa