大气海洋数据同化方法松弛逼近Kalman滤波入门课件

回顾逐步订正：背景场（猜值场、预报场、诊断场）+观测数据=优化场优点：简单，计算量小；观测数据较多时，效果还不错，不比最优插值逊色。缺点：只针对一个变量场没有动力与物理上的约束一个实例：由地面气压场诊断地面风场，然后利用地面风观测资料运用逐步订正法调诊断风场。

考虑惯性力（即平流项）、气压梯度力、科氏力、摩擦力，在局地直角平面建立一个简单的水平二维地面风场诊断模式：1980年4月4日20时的气压场诊断出地面/海面风场

订正前订正后

松弛逼近(Nudging,Newtonianrelaxation)其中，G为松弛系数，为观测向量。分析（格点）松弛逼近观测松弛逼近

最优松弛逼近局部效应全局效应待确定观测松弛逼近格点松弛逼近观测较少时观测较多时利用观测生成一个观测场

分析松弛法分析松弛法是对整个空间场的每一个格点都实施松弛调整。在实施分析松弛法之前，必须对观测（常常在背景场基础之上）进行客观分析从而得到观测格点空间场，它是与模式格点空间场一一对应的。利用整个客观分析场作为强迫场，其前提条件是客观分析场必须尽可能地反映真实天气形势。分析松弛法场适用于较大尺度的环境背景场的数据同化，常常用于常规气象观测资料的同化。如常规探空资料（每天两次，12小时一次）和常规地面观测资料（每3小时一次）。

格点松弛逼近(分析松弛逼近)

观测松弛法观测松弛法是对以观测为中心在其周围一定空间（即观测影响区域）之内的模式格点进行松弛调整，而在观测影响区域之外的区域不进行松弛调整。对同化非常规观测资料很有用，如探空加密观测资料，卫星资料和雷达资料等。

观测松弛逼近的空间影响范围

时间窗口观测松弛逼近的时间影响范围

一个实例：G是松弛系数，是观测质量系数（反映观测质量好坏和观测代表性）。是格点上的模式“预报”变量插值到观测点位置（与观测逼近）。

abca:诊断风场

b:逐步订正后c:观测松弛后

逐步订正与观测松弛逼近的区别在订正过程中随着订正次数的增加猜值风场一直无条件地向观测值逼近，既不考虑观测质量好坏，也不考虑物理上的动力约束。订正结果完全取决于观测值与订正次数的选取。观测误差不仅造成观测与猜值场之间的不协调，也造成不同观测之间的不协调。而且逐步订正是对风速两个水平分量分别独立进行，两个分量之间不仅由于没有物理上的动力约束而相互限制，且由于前面的不协调会彼此影响而进一步远离动力约束，因而最终有可能会导致调整场形势与猜值场大相径庭，虽然调整场在观测点与观测十分吻合。松弛逼近法不同，首先对观测进行粗略“评估”（即确定观测系数），然后让猜值风场在向观测值逼近的同时还要受限于物理模式的动力约束，大大降低那些不满足动力约束条件而与物理系统不相协调的观测资料在调整过程中对猜值风场的影响，让质量较高的观测发挥主要作用。

什么情况下运用逐步订正与观测松弛逼近？无动力约束方程，只能用逐步订正。反之，就可以用观测松弛逼近！松弛逼近的其它用途问题：嵌套网格计算时如何给嵌套区提供侧边界条件？finegrid边界条件的给定一个运用观测松弛逼近的简单例子

Kalman滤波入门Kalman滤波入门

20世纪60年代初期，Kalman（1960），Kalman和Bucy（1961）针对随机过程状态估计提出了Kalman滤波的思想。之后，Kalman滤波在信号处理、最优控制、航天等领域得到广泛应用。

20世纪60年代中期，Jones（1965）首次将Kalman滤波引入气象学。基本知识：dynamicsystem:x(t+1)=A*x(t)+u(k)measurementsystem:y(t)=M*x(t)+w(k)假设：

E(u)=0;E(u*u)=QE(w)=0;E(w*w)=R

高斯型白噪音（自相关）Kalman滤波示意图名词：

A状态转换距阵

M观测转换距阵状态(相）状态空间（相空间）为什么需要观测转换距阵？原因1：从观测点到模式点原因2：从观测变量到模式变量dynamicsystem:x(t+1)=A*x(t)+u(k)measurementsystem:y(t)=M*x(t)+w(k)Kalman滤波基本思想：

预报场初始场分析场=预报场+K（观测-模式解）最小方差估计P，K基本思路讲解：1)在t0时刻，给初始状态x(t0)一个估计xe(t0),误差协方差为PP(t0)=E[x(t0)-xe(t0)]**22)我们想得到x(t1)x(t1)=A*x(t0)+u(t0)xe(t1)=A*xe(t0)

P(t1)=E[x(t1)-xe(t1)]**2=A*P(t0)**2+Q

3)t1时刻观测系统y(t1)=M*x(t1)+w(t1)

在使用观测之前，ye(t1)=M*xe(t1)

使用观测之后，

newxe(t1)=xe(t1)+K[y(t1)-ye(1)]K:Kalmangain

于是，我们可以得到新P(t1),newP(t1)=E{x(t1)-newxe(t1)}**2=E{x(t1)-xe(t1)-K[y(t1)-ye(1)]}**2=E{[x(t1)-xe(t1)](1-k*M)-w*K}**2=P(t1)*(1-K*M)**2+R*K**2

得到，K=P*M/(P*M**2+R)

因此，newxe(t1)=A*xe(t1)+K*[y(t1)–ye(t1)]

然后

xe(t1)=newxe(t1)P(t1)=newP(t1)

重复以上过程往后走(t2,t3,t4,t5,………)！4)使newP(t1)最小，一个非常简单的例子：

x(t+1)=0.9*x(t)Q=100y(t+1)=y(t)R=10000初值：xe(t0)=1000,P(t0)=40000如果t1时刻有观测y(1)=1200,请估计x(t1)=?t0t1t2xe(t0)=1000xe(t1)=900P(t0)=40000K=0.7647y(1)