四维变分资料同化

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
四维变分资料同化
( 4DVAR, four-dimensional variational data assimilation method )
兰州大学大气科学学院 邱崇践
1,有关资料同化的基本知识
资料同化: 资料同化 在积分描写动力系统演变过程的数学模 预报模式)的同时,不断吸收观测资料, 式(预报模式)的同时,不断吸收观测资料,给出 系统状况的一个估计. 系统状况的一个估计.
(δ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
( 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) )
∫∫ ∫∫
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)
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) 向量空间内积: 向量空间内积:
b
∑
k =1
k k
代入目标函数可以直接求梯度
1 T J (α) = (M 1)αT ∧2 α + ∑[H( xB +V α) y]r R1[H( xB +V α) y]r 2 r
再
见
�
以 δ~为初值反向积分伴随模式到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) ) 伴随模式的解析形式只能作为理论推导用, 伴随模式的解析形式只能作为理论推导用,实际 问题是离散化的, 问题是离散化的,预报程序中还有些是不能写成 解析公式的.要保证相应的伴随模式严格成立, 解析公式的.要保证相应的伴随模式严格成立, 通常的作法是先根据原模式计算程序写出切线性 模式程序, 模式程序,再直接根据切线性模式程序一一对应 地写出伴随程序. 地写出伴随程序.一个天气预报模式的程序有上 万条语句,首先写出他的切线性模式程序, 万条语句,首先写出他的切线性模式程序,然后 根据切线性模式程序写出伴随程序, 根据切线性模式程序写出伴随程序,工作量是巨 大的. 大的. 按照一定的规则来写. 按照一定的规则来写.
(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*)
(需要梯度) 需要梯度)
最速下降法,共轭梯度法,拟牛顿法 最速下降法,共轭梯度法,拟牛顿法…
得到搜寻方向后成为一维寻优. 得到搜寻方向后成为一维寻优. 多项式逼近: 多项式逼近:
2,变分方法和伴随模式 , 变分方法是求泛函极值的有力工具. 变分方法是求泛函极值的有力工具. 是求泛函极值的有力工具 泛函的一个例子: 泛函的一个例子:
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的评论 的评论 好处:利用完整的模式方程约束,同化场动力上 好处:利用完整的模式方程约束, 协调,可以同时同化多时刻资料, 协调,可以同时同化多时刻资料,容易加 入其它约束.已经成功运用. 入其它约束.已经成功运用. 问题:不考虑模式误差; 问题:不考虑模式误差;
~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
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) )
示意图: 示意图:
如何求极小?下降算法 如何求极小?下降算法.
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的线性部分 阶 变分:
先看连续情况.反演初值的一个例子: 先看连续情况.反演初值的一个例子:目标泛函
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)
给出的最好估计是什么? 综合观测和背景场 给出的最好估计是什么? 背景场 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
目的:给出大气,海洋,陆面…状态的最好估计 状态的最好估计, 目的:给出大气,海洋,陆面 状态的最好估计, 为预报和分析研究提供必要的数据. 为预报和分析研究提供必要的数据.
为什么要用预报模式? 为什么要用预报模式? (1)观测不足 ) (2)观测有误差 ) 后果:变量间的不协调造成预报的振荡) (后果:变量间的不协调造成预报的振荡) 预报模式给我们提供什么? 预报模式给我们提供什么? (1)模式作出的预报为同化提供初猜场(背景场) )模式作出的预报为同化提供初猜场(背景场) (2)模式在不同点的变量之间以及各个变量之间建 ) 立了联系
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) 维变分方法: 极小化( ) 维变分方法