修正误差扰动的集合Kalman滤波方法在近地表土壤湿度同

1、田向军田向军 谢正辉谢正辉 师春香师春香中国科学院大气物理研究所中国科学院大气物理研究所ECCE Summer School for Advanced Study in Climate and Earth EnvironmentJuly 30-August 12, 2006, Beijing, China 提纲提纲 引言引言 集合集合Kalman滤波算法的实现滤波算法的实现 简单的土壤水模型简单的土壤水模型 理想验证试验理想验证试验 正在进行的相关工作正在进行的相关工作 一一、引言引言 土壤湿度在陆气相互作用的重要性土壤湿度在陆气相互作用的重要性 集合集合Kalman滤波方法在数据同化中的优劣

2、势滤波方法在数据同化中的优劣势二、集合二、集合Kalman滤波滤波 集合集合Kalman滤波是由滤波是由Evensen(1994b)所开发的所开发的一种数据同化方法一种数据同化方法.其详尽的说明请参阅其详尽的说明请参阅其相关其相关的文章的文章集合协方差矩阵集合协方差矩阵P 我们首先定义下面的集合状态矩阵我们首先定义下面的集合状态矩阵 定义定义 定义集合扰动矩阵定义集合扰动矩阵12(,)n NNAR A1NAA(1 )NAAAA I 由此由此,集合协方差矩阵集合协方差矩阵 便可以定义为便可以定义为n nePR()1TeAAPN观测扰动观测扰动 给定一个给定的观测向量给定一个给定的观测向量 (默认

3、的最优观默认的最优观测测). 是观测的个数是观测的个数. 定义定义 个具有扰动的观个具有扰动的观测向量如下测向量如下mdRmN,1,jjddjN 12(,)m NNDd ddR12( ,)m NNER 1TeEERN分析方程分析方程 利用集合斜方差矩阵可以将分析方程写成下面利用集合斜方差矩阵可以将分析方程写成下面的形式的形式 集合增量集合增量(innovation)向量定义为向量定义为1() ()aTTeeeAAPHHPHRDHADDHA 分析方程可以变为分析方程可以变为 1()aTTTTTAA AA H HAA HEED 引入矩阵引入矩阵 ,同时定义矩阵同时定义矩阵m NS HA Rm mC

4、 R(1)TCSSNR(1)TeeCSSNRTTSSEE1aTTAA AS S C D1(1 )TNA A IS C D1()TA IS C DAX1TXI S C D 集合集合Kalman滤波的一种平方根算滤波的一种平方根算法法 Evensen(2005)提出了一种新的平方根分析方法提出了一种新的平方根分析方法. 该方法能够避免观测引入的观测扰动对分析该方法能够避免观测引入的观测扰动对分析结果的影响结果的影响,另外该算法在计算分析的结构的时另外该算法在计算分析的结构的时候不需要引入任何额外的附加假设与近似候不需要引入任何额外的附加假设与近似,例如例如不需要假设观测扰动与集合状态标量扰动的无不

5、需要假设观测扰动与集合状态标量扰动的无关性关性,也不需要对观测误差斜方差矩阵求逆也不需要对观测误差斜方差矩阵求逆,其其算法简便有效算法简便有效. 这种新的算法是从传统的这种新的算法是从传统的Kalman滤波分析方滤波分析方程的斜方差矩阵的更新出发程的斜方差矩阵的更新出发,进行矩阵变化进行矩阵变化,得得到集合状态变量扰动的更新到集合状态变量扰动的更新 1()affTfTfPPP HHP HRHP1()aaTTTA AA I S C S A1()afTfAS CdH算法执行算法执行 1. 计算出计算出C C,然后对然后对C C进行特征值分解进行特征值分解 2.2.更新状态变量的集合平均值更新状态变

6、量的集合平均值 TZ ZC1()afTfAS ZdH按下列顺序进行矩阵乘法按下列顺序进行矩阵乘法 1()TfyZdH121yy32yZy43TyS y4afAy 3.3.计算矩阵计算矩阵 4.4.进行进行SVDSVD分解分解 5.5.求解状态变量的分析集合扰动求解状态变量的分析集合扰动 6.6. 然后再加上集合分析平均值然后再加上集合分析平均值 ,便得到状态便得到状态变量的分析扰动变量的分析扰动. 122TXZ S22 22TUVX222aTAAV Ia初始集合的生成初始集合的生成 如果初始样本的选择适当将有助于准确的表征如果初始样本的选择适当将有助于准确的表征状态初始量的误差特征状态初始量的

7、误差特征.但是但是,如果选择的样本如果选择的样本偏差不是特别大偏差不是特别大,随着时间的积分一般不会对同随着时间的积分一般不会对同化结果产生影响化结果产生影响.初始的样本一般采取这样的产初始的样本一般采取这样的产生方式生方式:首先对初始状态变量给出一个较为准首先对初始状态变量给出一个较为准确的估计值确的估计值,然后再在这个估计值的附近进行扰然后再在这个估计值的附近进行扰动动.产生扰动的方式有很多种。最简单的是产生产生扰动的方式有很多种。最简单的是产生随机数随机数(对标量模型对标量模型),随机曲线随机曲线(对对1维模型维模型),或或者是随机场者是随机场(对对2维或者是以上维数的模型维或者是以上维

8、数的模型),具具体模型的产生方法参加体模型的产生方法参加EnvensenEnvensen的文章的文章. 集合的积分集合的积分 对于一个离散非线性动力模型对于一个离散非线性动力模型1()kkkfW三三、用于验证模型的简单的土壤、用于验证模型的简单的土壤水模型水模型( )( )KDtzzz23( )bssKK2( )bssssbKD( , )sL t 四四、 理想验证试验理想验证试验简单验证试验简单验证试验05101520253035195189283377471565659753847941103511291223131714111505159916931787同化结果观测值模拟值真实值0510

9、15202530351815222936435057647178859299系列1系列2系列3系列40510152025303511835526986103120137154171188205222239256系列1系列2系列3系列405101520253015913172125293337414549系列1系列2系列3系列4051015202530147101316192225283134374043系列1系列2系列3系列405101520253035135791113151719212325系列1系列2系列3系列4 对上面简单的土壤水模型对上面简单的土壤水模型, 模型参数、初边值模型参数、

10、初边值、时间步长以及空间步长取下面的值、时间步长以及空间步长取下面的值: 0.54s63.2 10/sKms1800ts 0.1zm 2.0sm0.1/100/3600 /qm s00.54 0.4190.54 0.419r 利用该模型进行积分利用该模型进行积分,积分积分1800步作为模式的模拟值步作为模式的模拟值. 理想的观测值产生方式如下理想的观测值产生方式如下: 在上面的模型积分中在上面的模型积分中,在在每一个时间步随机产生一个均值为每一个时间步随机产生一个均值为0,标准差为标准差为0.01的随的随机扰动加到机扰动加到 上上, 也就是在积分的每一步都对土壤水也就是在积分的每一步都对土壤水

11、分含量的边界值产生一个随机扰动分含量的边界值产生一个随机扰动,由此同样积分由此同样积分1800步步,产生一套数据产生一套数据, 把这套数据作为观测数据把这套数据作为观测数据. 由于实际由于实际观测中观测中,我们往往只能得到陆地近地表的土壤湿度观测我们往往只能得到陆地近地表的土壤湿度观测值值. 因此因此,我们只选取每一步所产生的土壤湿度廓线靠我们只选取每一步所产生的土壤湿度廓线靠近地面的上面四层土壤的土壤湿度作为每一步的观测近地面的上面四层土壤的土壤湿度作为每一步的观测数据数据.(0)一般通用的集合一般通用的集合Kalman滤波方法滤波方法的应用的应用 利用上面描述的模型，以及上面构造的理想观利

12、用上面描述的模型，以及上面构造的理想观测数据进行同化试验测数据进行同化试验,其中模型的随机误差取为其中模型的随机误差取为均值为均值为0,标准差为标准差为0.03的随机向量的随机向量.下面给出的下面给出的是模拟的结果是模拟的结果.0200400600800100012001400160018000.220.240.260.280.30.320.340.360.38ObservationCommon methodSimulation第一层土壤观测模拟与同化的结果第一层土壤观测模拟与同化的结果0200400600800100012001400160018000.220.240.260.280.30.

13、320.340.36ObservationCommon methodSimulation第二层土壤观测模拟与同化的结果第二层土壤观测模拟与同化的结果0200400600800100012001400160018000.220.240.260.280.30.320.340.36ObservationCommon methodSimulation第三层土壤观测模拟与同化的结果第三层土壤观测模拟与同化的结果0200400600800100012001400160018000.220.240.260.280.30.320.34ObservationCommon methodSimulation第四层土

14、壤观测模拟与同化的结果第四层土壤观测模拟与同化的结果0200400600800100012001400160018000.220.230.240.250.260.270.280.290.3ObservationCommon methodSimulation第五层土壤观测模拟与同化的结果第五层土壤观测模拟与同化的结果0200400600800100012001400160018000.220.230.240.250.260.270.280.290.30.31ObservationCommon methodSimulation第六层土壤观测模拟与同化的结果第六层土壤观测模拟与同化的结果020040

15、0600800100012001400160018000.20.250.30.350.40.450.50.55ObservationCommon methodSimulation第三十层土壤观测模拟与同化的结果第三十层土壤观测模拟与同化的结果510152025300.250.30.350.40.450.5ObservationCommon methodSimulation一般的集合一般的集合Kalman滤波在第滤波在第600步的同化步的同化与观测土壤湿度廓线的对比与观测土壤湿度廓线的对比 到底是什么原因造成这么大的误差呢到底是什么原因造成这么大的误差呢?到底是到底是什么原因使得同化后的土壤湿度

16、廓线如此的奇什么原因使得同化后的土壤湿度廓线如此的奇怪怪? 这是一个问题这是一个问题!一些猜测一些猜测 问题的关键在于我们只有上面问题的关键在于我们只有上面4层的土壤湿度观测资料层的土壤湿度观测资料,而在生成集合对模式进行扰动的时候而在生成集合对模式进行扰动的时候,整个土壤整个土壤30层层(用于本文的这个简单的土壤水模型用于本文的这个简单的土壤水模型)都给予了都给予了0.03标准标准差的扰动差的扰动,而没有观测数据进来而没有观测数据进来,使得土壤廓线的曲线比使得土壤廓线的曲线比较怪异较怪异,在海洋模型中对于这些缺乏数据的点在海洋模型中对于这些缺乏数据的点,如果是一如果是一个平面上一般采用插值的

17、方法个平面上一般采用插值的方法,而我们这里显然不可以而我们这里显然不可以进行直接的插值进行直接的插值,当然有人采用根据某些特定的经验关当然有人采用根据某些特定的经验关系给出下面土壤湿度的系给出下面土壤湿度的“观测值观测值”. 而在这里既然没有而在这里既然没有观测我们就充分相信模式观测我们就充分相信模式,同时利用上面的观测信息去同时利用上面的观测信息去修正下面的廓线修正下面的廓线. 改进模型误差扰动的集合改进模型误差扰动的集合Kalman滤波方法的应用滤波方法的应用 一般的集合集合一般的集合集合Kalman滤波同化失败的原因发现，在滤波同化失败的原因发现，在该试验中该试验中(实际情况中也多是如此

18、实际情况中也多是如此),只有土壤表层的观只有土壤表层的观测值测值,也就是说只有在土壤表层有观测信息也就是说只有在土壤表层有观测信息,对于土壤其对于土壤其他层次的则没有观测信息他层次的则没有观测信息.而对于一般的集合而对于一般的集合Kalman滤波而言滤波而言,模型的模式误差一般取为均值为模型的模式误差一般取为均值为0,方差为方差为 的随机误差向量的随机误差向量(其维数与状态标量向量的维数一致其维数与状态标量向量的维数一致),也就是说对状态向量的每一个变量都进行随机扰动也就是说对状态向量的每一个变量都进行随机扰动,在在我们的这个数值验证中亦是如此我们的这个数值验证中亦是如此,因此对于深层次的土因

19、此对于深层次的土壤湿度而言壤湿度而言,较大随机误差的扰动同时缺乏必要的观测较大随机误差的扰动同时缺乏必要的观测信息信息,从而使得深层次的同化结果非常糟糕从而使得深层次的同化结果非常糟糕.既然缺乏充既然缺乏充分的土壤深层次观测信息分的土壤深层次观测信息, 我们只能充分相信模式的模拟结果我们只能充分相信模式的模拟结果,同时土壤表同时土壤表层的观测信息也必然通过模式的时间积分去影层的观测信息也必然通过模式的时间积分去影响深层次的土壤水分含量响深层次的土壤水分含量.鉴于此鉴于此,我们改变通我们改变通常使用的模式随机误差扰动的产生方式常使用的模式随机误差扰动的产生方式,只在有只在有观测信息的区域观测信息

20、的区域,比如前比如前m层层(m是有观测信息的是有观测信息的层数层数),产生随机误差扰动产生随机误差扰动,而对于而对于m+1到到n层认层认为误差为为误差为0.下面是改变了模式随机误差扰动后下面是改变了模式随机误差扰动后的同化结果分析。的同化结果分析。下面是改进误差扰动方式的下面是改进误差扰动方式的模拟的结果模拟的结果0200400600800100012001400160018000.220.240.260.280.30.320.340.360.38ObservationNew methodOld methodSimulation第一层土壤观测模拟与同化的结第一层土壤观测模拟与同化的结果果020

21、0400600800100012001400160018000.220.240.260.280.30.320.340.36ObservationNew methodOld methodSimulation第二层土壤观测模拟与同化的结果第二层土壤观测模拟与同化的结果0200400600800100012001400160018000.220.240.260.280.30.320.340.36ObservationNew methodOld methodSimulation第三层土壤观测模拟与同化的结果第三层土壤观测模拟与同化的结果0200400600800100012001400160018000.220.240.260.280.30.320.34ObservationNew methodOld methodSimulation第四层土壤观测模拟与同化的结果第四层土壤观测模拟与同化的结果510152025300.250.30.350.40.450.5ObservationNew methodOld meth