修正誤差擾動的集合Kalman濾波方法在近地表土壤濕度同

1、田向軍田向軍 謝正輝謝正輝 師春香師春香中國科學院大氣物理研究所中國科學院大氣物理研究所ECCE Summer School for Advanced Study in Climate and Earth EnvironmentJuly 30-August 12, 2006, Beijing, China 提綱提綱 引言引言 集合集合Kalman濾波算法的實現(xiàn)濾波算法的實現(xiàn) 簡單的土壤水模型簡單的土壤水模型 理想驗證試驗理想驗證試驗 正在進行的相關(guān)工作正在進行的相關(guān)工作 一一、引言引言 土壤濕度在陸氣相互作用的重要性土壤濕度在陸氣相互作用的重要性 集合集合Kalman濾波方法在數(shù)據(jù)同化中的優(yōu)劣

2、勢濾波方法在數(shù)據(jù)同化中的優(yōu)劣勢二、集合二、集合Kalman濾波濾波 集合集合Kalman濾波是由濾波是由Evensen(1994b)所開發(fā)的所開發(fā)的一種數(shù)據(jù)同化方法一種數(shù)據(jù)同化方法.其詳盡的說明請參閱其詳盡的說明請參閱其相關(guān)其相關(guān)的文章的文章集合協(xié)方差矩陣集合協(xié)方差矩陣P 我們首先定義下面的集合狀態(tài)矩陣我們首先定義下面的集合狀態(tài)矩陣 定義定義 定義集合擾動矩陣定義集合擾動矩陣12(,)n NNAR A1NAA(1 )NAAAA I 由此由此,集合協(xié)方差矩陣集合協(xié)方差矩陣 便可以定義為便可以定義為n nePR()1TeAAPN觀測擾動觀測擾動 給定一個給定的觀測向量給定一個給定的觀測向量 (默認

3、的最優(yōu)觀默認的最優(yōu)觀測測). 是觀測的個數(shù)是觀測的個數(shù). 定義定義 個具有擾動的觀個具有擾動的觀測向量如下測向量如下mdRmN,1,jjddjN 12(,)m NNDd ddR12( ,)m NNER 1TeEERN分析方程分析方程 利用集合斜方差矩陣可以將分析方程寫成下面利用集合斜方差矩陣可以將分析方程寫成下面的形式的形式 集合增量集合增量(innovation)向量定義為向量定義為1() ()aTTeeeAAPHHPHRDHADDHA 分析方程可以變?yōu)榉治龇匠炭梢宰優(yōu)?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)提出了一種新的平方根分析方法提出了一種新的平方根分析方法. 該方法能夠避免觀測引入的觀測擾動對分析該方法能夠避免觀測引入的觀測擾動對分析結(jié)果的影響結(jié)果的影響,另外該算法在計算分析的結(jié)構(gòu)的時另外該算法在計算分析的結(jié)構(gòu)的時候不需要引入任何額外的附加假設與近似候不需要引入任何額外的附加假設與近似,例如例如不需要假設觀測擾動與集合狀態(tài)標量擾動的無不

5、需要假設觀測擾動與集合狀態(tài)標量擾動的無關(guān)性關(guān)性,也不需要對觀測誤差斜方差矩陣求逆也不需要對觀測誤差斜方差矩陣求逆,其其算法簡便有效算法簡便有效. 這種新的算法是從傳統(tǒng)的這種新的算法是從傳統(tǒng)的Kalman濾波分析方濾波分析方程的斜方差矩陣的更新出發(fā)程的斜方差矩陣的更新出發(fā),進行矩陣變化進行矩陣變化,得得到集合狀態(tài)變量擾動的更新到集合狀態(tài)變量擾動的更新 1()affTfTfPPP HHP HRHP1()aaTTTA AA I S C S A1()afTfAS CdH算法執(zhí)行算法執(zhí)行 1. 計算出計算出C C,然后對然后對C C進行特征值分解進行特征值分解 2.2.更新狀態(tài)變量的集合平均值更新狀態(tài)變

6、量的集合平均值 TZ ZC1()afTfAS ZdH按下列順序進行矩陣乘法按下列順序進行矩陣乘法 1()TfyZdH121yy32yZy43TyS y4afAy 3.3.計算矩陣計算矩陣 4.4.進行進行SVDSVD分解分解 5.5.求解狀態(tài)變量的分析集合擾動求解狀態(tài)變量的分析集合擾動 6.6. 然后再加上集合分析平均值然后再加上集合分析平均值 ,便得到狀態(tài)便得到狀態(tài)變量的分析擾動變量的分析擾動. 122TXZ S22 22TUVX222aTAAV Ia初始集合的生成初始集合的生成 如果初始樣本的選擇適當將有助于準確的表征如果初始樣本的選擇適當將有助于準確的表征狀態(tài)初始量的誤差特征狀態(tài)初始量的

7、誤差特征.但是但是,如果選擇的樣本如果選擇的樣本偏差不是特別大偏差不是特別大,隨著時間的積分一般不會對同隨著時間的積分一般不會對同化結(jié)果產(chǎn)生影響化結(jié)果產(chǎn)生影響.初始的樣本一般采取這樣的產(chǎn)初始的樣本一般采取這樣的產(chǎn)生方式生方式:首先對初始狀態(tài)變量給出一個較為準首先對初始狀態(tài)變量給出一個較為準確的估計值確的估計值,然后再在這個估計值的附近進行擾然后再在這個估計值的附近進行擾動動.產(chǎn)生擾動的方式有很多種。最簡單的是產(chǎn)生產(chǎn)生擾動的方式有很多種。最簡單的是產(chǎn)生隨機數(shù)隨機數(shù)(對標量模型對標量模型),隨機曲線隨機曲線(對對1維模型維模型),或或者是隨機場者是隨機場(對對2維或者是以上維數(shù)的模型維或者是以上維

8、數(shù)的模型),具具體模型的產(chǎn)生方法參加體模型的產(chǎn)生方法參加EnvensenEnvensen的文章的文章. 集合的積分集合的積分 對于一個離散非線性動力模型對于一個離散非線性動力模型1()kkkfW三三、用于驗證模型的簡單的土壤、用于驗證模型的簡單的土壤水模型水模型( )( )KDtzzz23( )bssKK2( )bssssbKD( , )sL t 四四、 理想驗證試驗理想驗證試驗簡單驗證試驗簡單驗證試驗05101520253035195189283377471565659753847941103511291223131714111505159916931787同化結(jié)果觀測值模擬值真實值0510

9、15202530351815222936435057647178859299系列1系列2系列3系列40510152025303511835526986103120137154171188205222239256系列1系列2系列3系列405101520253015913172125293337414549系列1系列2系列3系列4051015202530147101316192225283134374043系列1系列2系列3系列405101520253035135791113151719212325系列1系列2系列3系列4 對上面簡單的土壤水模型對上面簡單的土壤水模型, 模型參數(shù)、初邊值模型參數(shù)、

10、初邊值、時間步長以及空間步長取下面的值、時間步長以及空間步長取下面的值: 0.54s63.2 10/sKms1800ts 0.1zm 2.0sm0.1/100/3600 /qm s00.54 0.4190.54 0.419r 利用該模型進行積分利用該模型進行積分,積分積分1800步作為模式的模擬值步作為模式的模擬值. 理想的觀測值產(chǎn)生方式如下理想的觀測值產(chǎn)生方式如下: 在上面的模型積分中在上面的模型積分中,在在每一個時間步隨機產(chǎn)生一個均值為每一個時間步隨機產(chǎn)生一個均值為0,標準差為標準差為0.01的隨的隨機擾動加到機擾動加到 上上, 也就是在積分的每一步都對土壤水也就是在積分的每一步都對土壤水

11、分含量的邊界值產(chǎn)生一個隨機擾動分含量的邊界值產(chǎn)生一個隨機擾動,由此同樣積分由此同樣積分1800步步,產(chǎn)生一套數(shù)據(jù)產(chǎn)生一套數(shù)據(jù), 把這套數(shù)據(jù)作為觀測數(shù)據(jù)把這套數(shù)據(jù)作為觀測數(shù)據(jù). 由于實際由于實際觀測中觀測中,我們往往只能得到陸地近地表的土壤濕度觀測我們往往只能得到陸地近地表的土壤濕度觀測值值. 因此因此,我們只選取每一步所產(chǎn)生的土壤濕度廓線靠我們只選取每一步所產(chǎn)生的土壤濕度廓線靠近地面的上面四層土壤的土壤濕度作為每一步的觀測近地面的上面四層土壤的土壤濕度作為每一步的觀測數(shù)據(jù)數(shù)據(jù).(0)一般通用的集合一般通用的集合Kalman濾波方法濾波方法的應用的應用 利用上面描述的模型，以及上面構(gòu)造的理想觀利

12、用上面描述的模型，以及上面構(gòu)造的理想觀測數(shù)據(jù)進行同化試驗測數(shù)據(jù)進行同化試驗,其中模型的隨機誤差取為其中模型的隨機誤差取為均值為均值為0,標準差為標準差為0.03的隨機向量的隨機向量.下面給出的下面給出的是模擬的結(jié)果是模擬的結(jié)果.0200400600800100012001400160018000.220.240.260.280.30.320.340.360.38ObservationCommon methodSimulation第一層土壤觀測模擬與同化的結(jié)果第一層土壤觀測模擬與同化的結(jié)果0200400600800100012001400160018000.220.240.260.280.30.

13、320.340.36ObservationCommon methodSimulation第二層土壤觀測模擬與同化的結(jié)果第二層土壤觀測模擬與同化的結(jié)果0200400600800100012001400160018000.220.240.260.280.30.320.340.36ObservationCommon methodSimulation第三層土壤觀測模擬與同化的結(jié)果第三層土壤觀測模擬與同化的結(jié)果0200400600800100012001400160018000.220.240.260.280.30.320.34ObservationCommon methodSimulation第四層土

14、壤觀測模擬與同化的結(jié)果第四層土壤觀測模擬與同化的結(jié)果0200400600800100012001400160018000.220.230.240.250.260.270.280.290.3ObservationCommon methodSimulation第五層土壤觀測模擬與同化的結(jié)果第五層土壤觀測模擬與同化的結(jié)果0200400600800100012001400160018000.220.230.240.250.260.270.280.290.30.31ObservationCommon methodSimulation第六層土壤觀測模擬與同化的結(jié)果第六層土壤觀測模擬與同化的結(jié)果020040

15、0600800100012001400160018000.20.250.30.350.40.450.50.55ObservationCommon methodSimulation第三十層土壤觀測模擬與同化的結(jié)果第三十層土壤觀測模擬與同化的結(jié)果510152025300.250.30.350.40.450.5ObservationCommon methodSimulation一般的集合一般的集合Kalman濾波在第濾波在第600步的同化步的同化與觀測土壤濕度廓線的對比與觀測土壤濕度廓線的對比 到底是什么原因造成這么大的誤差呢到底是什么原因造成這么大的誤差呢?到底是到底是什么原因使得同化后的土壤濕度

16、廓線如此的奇什么原因使得同化后的土壤濕度廓線如此的奇怪怪? 這是一個問題這是一個問題!一些猜測一些猜測 問題的關(guān)鍵在于我們只有上面問題的關(guān)鍵在于我們只有上面4層的土壤濕度觀測資料層的土壤濕度觀測資料,而在生成集合對模式進行擾動的時候而在生成集合對模式進行擾動的時候,整個土壤整個土壤30層層(用于本文的這個簡單的土壤水模型用于本文的這個簡單的土壤水模型)都給予了都給予了0.03標準標準差的擾動差的擾動,而沒有觀測數(shù)據(jù)進來而沒有觀測數(shù)據(jù)進來,使得土壤廓線的曲線比使得土壤廓線的曲線比較怪異較怪異,在海洋模型中對于這些缺乏數(shù)據(jù)的點在海洋模型中對于這些缺乏數(shù)據(jù)的點,如果是一如果是一個平面上一般采用插值的

17、方法個平面上一般采用插值的方法,而我們這里顯然不可以而我們這里顯然不可以進行直接的插值進行直接的插值,當然有人采用根據(jù)某些特定的經(jīng)驗關(guān)當然有人采用根據(jù)某些特定的經(jīng)驗關(guān)系給出下面土壤濕度的系給出下面土壤濕度的“觀測值觀測值”. 而在這里既然沒有而在這里既然沒有觀測我們就充分相信模式觀測我們就充分相信模式,同時利用上面的觀測信息去同時利用上面的觀測信息去修正下面的廓線修正下面的廓線. 改進模型誤差擾動的集合改進模型誤差擾動的集合Kalman濾波方法的應用濾波方法的應用 一般的集合集合一般的集合集合Kalman濾波同化失敗的原因發(fā)現(xiàn)，在濾波同化失敗的原因發(fā)現(xiàn)，在該試驗中該試驗中(實際情況中也多是如此

18、實際情況中也多是如此),只有土壤表層的觀只有土壤表層的觀測值測值,也就是說只有在土壤表層有觀測信息也就是說只有在土壤表層有觀測信息,對于土壤其對于土壤其他層次的則沒有觀測信息他層次的則沒有觀測信息.而對于一般的集合而對于一般的集合Kalman濾波而言濾波而言,模型的模式誤差一般取為均值為模型的模式誤差一般取為均值為0,方差為方差為 的隨機誤差向量的隨機誤差向量(其維數(shù)與狀態(tài)標量向量的維數(shù)一致其維數(shù)與狀態(tài)標量向量的維數(shù)一致),也就是說對狀態(tài)向量的每一個變量都進行隨機擾動也就是說對狀態(tài)向量的每一個變量都進行隨機擾動,在在我們的這個數(shù)值驗證中亦是如此我們的這個數(shù)值驗證中亦是如此,因此對于深層次的土因

19、此對于深層次的土壤濕度而言壤濕度而言,較大隨機誤差的擾動同時缺乏必要的觀測較大隨機誤差的擾動同時缺乏必要的觀測信息信息,從而使得深層次的同化結(jié)果非常糟糕從而使得深層次的同化結(jié)果非常糟糕.既然缺乏充既然缺乏充分的土壤深層次觀測信息分的土壤深層次觀測信息, 我們只能充分相信模式的模擬結(jié)果我們只能充分相信模式的模擬結(jié)果,同時土壤表同時土壤表層的觀測信息也必然通過模式的時間積分去影層的觀測信息也必然通過模式的時間積分去影響深層次的土壤水分含量響深層次的土壤水分含量.鑒于此鑒于此,我們改變通我們改變通常使用的模式隨機誤差擾動的產(chǎn)生方式常使用的模式隨機誤差擾動的產(chǎn)生方式,只在有只在有觀測信息的區(qū)域觀測信息

20、的區(qū)域,比如前比如前m層層(m是有觀測信息的是有觀測信息的層數(shù)層數(shù)),產(chǎn)生隨機誤差擾動產(chǎn)生隨機誤差擾動,而對于而對于m+1到到n層認層認為誤差為為誤差為0.下面是改變了模式隨機誤差擾動后下面是改變了模式隨機誤差擾動后的同化結(jié)果分析。的同化結(jié)果分析。下面是改進誤差擾動方式的下面是改進誤差擾動方式的模擬的結(jié)果模擬的結(jié)果0200400600800100012001400160018000.220.240.260.280.30.320.340.360.38ObservationNew methodOld methodSimulation第一層土壤觀測模擬與同化的結(jié)第一層土壤觀測模擬與同化的結(jié)果果020

21、0400600800100012001400160018000.220.240.260.280.30.320.340.36ObservationNew methodOld methodSimulation第二層土壤觀測模擬與同化的結(jié)果第二層土壤觀測模擬與同化的結(jié)果0200400600800100012001400160018000.220.240.260.280.30.320.340.36ObservationNew methodOld methodSimulation第三層土壤觀測模擬與同化的結(jié)果第三層土壤觀測模擬與同化的結(jié)果0200400600800100012001400160018000.220.240.260.280.30.320.34ObservationNew methodOld methodSimulation第四層土壤觀測模擬與同化的結(jié)果第四層土壤觀測模擬與同化的結(jié)果510152025300.250.30.350.40.450.5ObservationNew methodOld meth