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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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