擺玉龍，尤元紅，邵　宇，張轉(zhuǎn)花

(西北師范大學(xué)物理與電子工程學(xué)院，甘肅蘭州　730070)

Lorenz混沌系統(tǒng)的數(shù)據(jù)同化方法研究

擺玉龍，尤元紅，邵宇，張轉(zhuǎn)花

(西北師范大學(xué)物理與電子工程學(xué)院，甘肅蘭州730070)

摘要：利用數(shù)據(jù)同化方法研究了Lorenz混沌系統(tǒng)中的非線性問題.提出了一種基于Cholesky分解的降秩平方根濾波算法，可以改善數(shù)據(jù)同化中的濾波發(fā)散現(xiàn)象.通過對(duì)卡爾曼濾波誤差協(xié)方差矩陣進(jìn)行Cholesky分解，降低協(xié)方差矩陣的計(jì)算量和存儲(chǔ)量，提高濾波的收斂速度.在Lorenz混沌系統(tǒng)上研究了降秩平方根濾波的性能，通過敏感性分析試驗(yàn)，討論了降秩平方根濾波的穩(wěn)定性，驗(yàn)證了算法的有效性，比較了降秩平方根濾波與集合卡爾曼濾波的同化性能.結(jié)果表明，在Lorenz混沌系統(tǒng)的短期預(yù)報(bào)實(shí)驗(yàn)中，降秩平方根濾波的同化性能優(yōu)于集合卡爾曼濾波.

關(guān)鍵詞：Cholesky分解；Lorenz混沌系統(tǒng)；降秩平方根濾波

收稿日期：2015-05-15;修改稿收到日期:2015-09-20

E-mail:yulongbai@gmail.com

基金項(xiàng)目:國(guó)家自然科學(xué)基金資助項(xiàng)目(41461078)

作者簡(jiǎn)介：擺玉龍(1973—)，男，甘肅會(huì)寧人，教授，碩士研究生導(dǎo)師.主要研究方向?yàn)閿?shù)據(jù)同化和參數(shù)估計(jì).

中圖分類號(hào)：TP 309.7

文獻(xiàn)標(biāo)志碼：標(biāo)志碼：A

文章編號(hào)：章編號(hào)：1001-988Ⅹ(2015)06-0044-06

Abstract:Data assimilation methods are used to study the nonlinear problems in Lorenz chaotic system.A reduced rank square root filter algorithm(RRSQRT) based on Cholesky decomposition is proposed,which can alleviate the filtering divergence phenomenon in data assimilation.Using Kalman filtering error covariance matrix Cholesky decomposition,it can reduce the amount of calculation and storage of covariance matrix and improve the convergent speed of filtering.With the Lorenz chaotic system,sensitivity analysis tests are conducted to research the performance of singular square root filter,the stability of the singular square root filter is discussed and the effectiveness of the algorithm is verified.The assimilation performance between RRSQRT and ensemble Kalman filter(EnKF) BFQis compared.The results show that the RRSQRT is better than EnKF in short time forecast experiment with Lorenz chaotic system.

Data assimilation methods research for Lorenz chaotic systems

BAI Yu-long，YOU Yuan-hong，SHAO Yu，ZHANG Zhuan-hua

(College of Physics and Electronic Engineering，Northwest Normal University，Lanzhou 730070，Gansu，China)

Key words:Cholesky decomposition;Lorenz chaotic systems;reduced rank square root filter

Lorenz以著名的Lorenz方程形式給出了混沌系統(tǒng)的實(shí)例[1].在耗散系統(tǒng)中，不同的初始條件下，確定性方程的解出現(xiàn)了混沌現(xiàn)象.混沌系統(tǒng)的研究目前包括混沌控制、混沌保密通信和混沌同步研究等[2-3].在混沌控制研究中，大多數(shù)的控制方法都是在已知參數(shù)的情形下給出的.如果系統(tǒng)的參數(shù)未知或者不確定，將會(huì)使控制方法失效.數(shù)據(jù)同化(也稱為資料同化)是源自大氣和海洋等領(lǐng)域的研究方法，其核心思想是在模式狀態(tài)預(yù)報(bào)中引入觀測(cè)數(shù)據(jù)，再根據(jù)觀測(cè)數(shù)據(jù)對(duì)模式狀態(tài)進(jìn)行更新，提高預(yù)報(bào)的效率和準(zhǔn)確性[4].卡爾曼濾波(Kalman Filter)方法是一種典型的順序數(shù)據(jù)同化方法[5].傳統(tǒng)卡爾曼濾波需要對(duì)其模型算子和觀測(cè)算子進(jìn)行線性化，很多研究者提出了卡爾曼濾波的擴(kuò)展形式，如擴(kuò)展卡爾曼濾波(EKF)、集合卡爾曼濾波(EnKF)[6]及確定性集合卡爾曼濾波(DEnKF)[7]等.針對(duì)EnKF分析過程中低估預(yù)報(bào)誤差協(xié)方差導(dǎo)致同化不能有效吸收觀測(cè)信息，造成濾波發(fā)散的問題，Whitaker[8]在2002年提出了無(wú)需擾動(dòng)觀測(cè)的集合平方根濾波(EnSRF).李昊睿[9]等提出了一種集合平方根濾波和四維變分(4DVAR)的混合方法來(lái)反演土壤濕度廓線，結(jié)果表明，混合方法反演的分析時(shí)刻土壤濕度廓線都優(yōu)于EnSRF和4DVAR的結(jié)果.王世璋等[10]在集合平方根濾波的基礎(chǔ)上提出了迭代集合平方根濾波，并成功應(yīng)用于風(fēng)暴尺度資料同化中.然而，利用數(shù)據(jù)同化方法研究混沌控制還較為少見.曹小群等[11]提出了基于MCMC方法的Lorenz混沌系統(tǒng)的參數(shù)估計(jì)，導(dǎo)出未知參數(shù)分布規(guī)律的后驗(yàn)概率密度函數(shù)；接著采用自適應(yīng)Metropolis算法構(gòu)造了Markov鏈；然后截取收斂的鏈序列，計(jì)算混沌系統(tǒng)參數(shù)的估計(jì)值.

文中提出一種新的數(shù)據(jù)同化方法，從卡爾曼濾波的協(xié)方差矩陣入手，在估計(jì)協(xié)方差矩陣的過程中引入Cholesky分解方法，將協(xié)方差矩陣分解為其平方根的乘積，并將這種變形的卡爾曼濾波結(jié)合Lorenz混沌系統(tǒng)進(jìn)行同化試驗(yàn).

1理論方法

1.1卡爾曼濾波

Kalman濾波作為一種順序數(shù)據(jù)同化濾波方法，其主要步驟為：首先對(duì)模式狀態(tài)進(jìn)行預(yù)報(bào)，接著引入觀測(cè)數(shù)據(jù)，然后利用觀測(cè)數(shù)據(jù)對(duì)模式狀態(tài)和誤差協(xié)方差進(jìn)行重新分析，得到變量的分析值.Kalman濾波的基本計(jì)算為：

1)預(yù)報(bào).利用i-1時(shí)刻狀態(tài)的分析值和誤差協(xié)方差的分析值通過模型預(yù)報(bào)，得到i時(shí)刻的預(yù)報(bào)值.

2)計(jì)算增益矩陣.利用i時(shí)刻的預(yù)報(bào)誤差協(xié)方差矩陣計(jì)算增益矩陣.

其中H為觀測(cè)矩陣；HT為觀測(cè)矩陣的轉(zhuǎn)置.

3)分析.引入觀測(cè)，對(duì)狀態(tài)量和誤差協(xié)方差進(jìn)行分析.

1.2降秩平方根濾波

降秩平方根濾波的基本思想是對(duì)狀態(tài)估計(jì)的協(xié)方差矩陣P進(jìn)行Cholesky分解，使協(xié)方差矩陣P=LLT，其中矩陣L具有q個(gè)特征值，且其對(duì)應(yīng)的特征向量為li,i=1,…,q,算法的執(zhí)行過程為：

1)隨機(jī)產(chǎn)生狀態(tài)變量.以任意的隨機(jī)數(shù)作為初始樣本，并以這組樣本對(duì)數(shù)據(jù)同化系統(tǒng)進(jìn)行初始化，同時(shí)導(dǎo)入觀測(cè)數(shù)據(jù)對(duì)樣本進(jìn)行實(shí)時(shí)更新.

其中，xa(0)為初始狀態(tài)變量分析值；La(0)為初始特征向量矩陣.

2)模型預(yù)報(bào).通過求解樣本狀態(tài)的控制方程得到狀態(tài)向量在第k個(gè)同化步的預(yù)報(bào)值，鑒于特征向量不是模型的狀態(tài)變量模型，不能直接利用模型對(duì)其進(jìn)行預(yù)報(bào)，現(xiàn)按(8)式將特征量投影到模型狀態(tài)變量中，再用模型進(jìn)行預(yù)報(bào)，系統(tǒng)的狀態(tài)方程為：

3)分析.

上式采用Cholesky分解的逆向思想，將模型預(yù)報(bào)獲得的平方根矩陣Lf(k)與其轉(zhuǎn)置矩陣Lf(k)T作積獲得協(xié)方差矩陣Pf(k).

圖1　降秩平方根濾波算法流程圖

2數(shù)值試驗(yàn)

2.1Lorenz混沌系統(tǒng)

Lorenz混沌系統(tǒng)最初用于描述容器底部加熱，容器內(nèi)液體的運(yùn)動(dòng)情況.容器底部液體越來(lái)越熱逐漸上升，形成對(duì)流，當(dāng)容器受熱到一定程度保持不變時(shí)，對(duì)流會(huì)以不規(guī)則湍流的方式運(yùn)動(dòng).這個(gè)系統(tǒng)經(jīng)過傅里葉分解、簡(jiǎn)化后，得到一個(gè)三維非線性方程組.Lorenz系統(tǒng)模型的數(shù)學(xué)表達(dá)式為：

(16)

其中,x為描述對(duì)流強(qiáng)度的量;y和z分別描述對(duì)流圈的水平和垂直方向溫度梯度[12]；σ和r分別與流體的Prandtl數(shù)和Rayleigh數(shù)成比例;b為與空間相關(guān)的常數(shù).當(dāng)系統(tǒng)參數(shù)σ=10,b=8/3,r=28時(shí)，系統(tǒng)表現(xiàn)出混沌現(xiàn)象，通常采用四階龍格庫(kù)塔算法求解系統(tǒng)微分方程.試驗(yàn)以Lorenz混沌系統(tǒng)為預(yù)報(bào)模型，以降秩平方根濾波為同化算法，對(duì)混沌系統(tǒng)的狀態(tài)變量進(jìn)行實(shí)時(shí)估計(jì).試驗(yàn)中采用均方根誤差RMSE評(píng)價(jià)模型狀態(tài)變量x估計(jì)值xe與真實(shí)值xt之間的差異，RMSE值越小，表明預(yù)報(bào)值與真實(shí)值越接近，同化效果越好.

其中N為模型狀態(tài)向量矩陣的大小.

2.2降秩平方根濾波基本同化試驗(yàn)

試驗(yàn)中取積分步長(zhǎng)h=0.025，設(shè)定積分區(qū)間為[0,10]，降秩平方根濾波協(xié)方差特征根數(shù)L=30，觀測(cè)間隔為1，觀測(cè)方差為0.5.僅取部分狀態(tài)變量值作為觀測(cè)量，分別抽取0，10 h，20 h，…,400 h時(shí)刻的狀態(tài)量作為n=40個(gè)觀測(cè)數(shù)據(jù)，在所選取的40個(gè)觀測(cè)數(shù)據(jù)上添加高斯隨機(jī)數(shù)作為觀測(cè)噪聲，即服從N(0,1).設(shè)定系統(tǒng)初始狀態(tài)值(x0,y0,z0)為(1.5088，-1.5312，25.4609)，圖2給出了混沌系統(tǒng)第400步時(shí)的基本同化結(jié)果.可見真值與估計(jì)值的大致走向相同，說明通過降秩平方根濾波預(yù)報(bào)的混沌系統(tǒng)動(dòng)力特性與混沌系統(tǒng)真實(shí)的動(dòng)力特性基本相同，表明降秩平方根濾波用于混沌系統(tǒng)的狀態(tài)預(yù)報(bào)是可行的.圖中實(shí)線與虛線距離越近，則預(yù)報(bào)值與真值之間的誤差越小，即同化效果越好.

圖2　第400步同化效果圖

2.3觀測(cè)方差對(duì)同化效果影響

為研究觀測(cè)誤差的變化對(duì)同化效果的影響，將試驗(yàn)中的觀測(cè)方差設(shè)置在[0.1,1.0]之間變化，對(duì)比不同觀測(cè)方差情況下的同化效果，分析觀測(cè)方差增大對(duì)降秩平方根濾波同化效果的影響.圖3給出了觀測(cè)噪聲在0.1～1.0之間變化時(shí)同化系統(tǒng)均方根誤差的變化情況.

從圖3可以看出，不同觀測(cè)方差條件下的RMSE也不同，隨著觀測(cè)方差的增大，狀態(tài)變量的RMSE值隨之增加.觀測(cè)方差為0.1時(shí)，RMSE值僅為0.12，觀測(cè)方差為0.2時(shí)，RMSE值為0.41，曲線的走向表明觀測(cè)方差越大，同化性能越差.這與應(yīng)用較為廣泛的集合卡爾曼濾波(EnKF)方法結(jié)論相一致.

圖3　觀測(cè)方差對(duì)同化效果的影響

2.4觀測(cè)窗口長(zhǎng)度對(duì)同化效果影響

為研究觀測(cè)窗口長(zhǎng)度對(duì)同化效果的影響，保持上述試驗(yàn)設(shè)置不變而僅改變同化窗口長(zhǎng)度，使得同化窗口長(zhǎng)度W分別取5,15,25，對(duì)比分析同化窗口長(zhǎng)度的改變對(duì)同化效果的影響.圖4給出了系統(tǒng)運(yùn)行400步時(shí)，同化窗口隨狀態(tài)變量x的RMSE值變化情況.不同的同化窗口長(zhǎng)度使得狀態(tài)變量x的RMSE值變化情況不盡相同.

1)觀測(cè)窗口W=5時(shí)，系統(tǒng)在運(yùn)行至300步之前，狀態(tài)變量x的均方根誤差保持在1以下，300步以后的RMSE值與300步之前相比雖有增大趨勢(shì)，但RMSE值總體較小，同化效果比較穩(wěn)定.

2)觀測(cè)窗口W=15時(shí)，系統(tǒng)運(yùn)行至130步后，誤差明顯增大，同化過程中的RMSE值與W=5時(shí)的RMSE值相比，有明顯變大趨勢(shì).

3)觀測(cè)窗口W=25時(shí)，系統(tǒng)運(yùn)行至130步后，濾波出現(xiàn)短暫發(fā)散現(xiàn)象，同化過程中的RMSE值與前兩個(gè)同化窗口長(zhǎng)度相比，W=25時(shí)的RMSE明顯要大很多.

圖4　觀測(cè)窗口對(duì)同化效果影響

實(shí)驗(yàn)表明，不同的同化窗口長(zhǎng)度對(duì)同化效果的影響不盡相同，同化窗口長(zhǎng)度越小，系統(tǒng)同化效果越好.三種同化效果的差異表明，同化系統(tǒng)中及時(shí)利用觀測(cè)數(shù)據(jù)更新模型的背景場(chǎng)，實(shí)時(shí)調(diào)整模型運(yùn)行軌跡，對(duì)提高同化效果有很大的幫助.

2.5加密觀測(cè)

考慮到濾波在局部出現(xiàn)較明顯的發(fā)散現(xiàn)象，文中嘗試采用加密觀測(cè)的方式對(duì)系統(tǒng)局部狀態(tài)預(yù)報(bào)進(jìn)行改善.從圖4可以看出，同化窗口為5時(shí)，系統(tǒng)在150步左右和300～400步的預(yù)報(bào)誤差較大.針對(duì)這種現(xiàn)象，試驗(yàn)在0～120步同化窗口長(zhǎng)度設(shè)定為5，120步以后采用觀測(cè)加密的措施，設(shè)定同化窗口為1，圖5給出了加密觀測(cè)前后的比較圖.

圖5　采用加密觀測(cè)前后比較

從圖5(a)可以看出，121～400步的預(yù)報(bào)值曲線在局部地方與真值的差值較大，與圖4中同化窗口W=5時(shí)反映的情況一致.圖5(b)表示在121～400步采用加密觀測(cè)手段后的效果圖，通過對(duì)比發(fā)現(xiàn)，預(yù)報(bào)值的曲線更接近于真值曲線，即表明加密觀測(cè)明顯改善了預(yù)報(bào)效果.

2.6與集合卡爾曼濾波(EnKF)同化性能比較

為進(jìn)一步研究降秩平方根濾波的同化性能，將降秩平方根濾波的性能與集合卡爾曼濾波的性能進(jìn)行比較.試驗(yàn)中采用并行計(jì)算技術(shù)，在Lorenz混沌系統(tǒng)上同時(shí)用兩種濾波方法對(duì)系統(tǒng)狀態(tài)進(jìn)行預(yù)報(bào)，降秩平方根濾波誤差協(xié)方差的特征根數(shù)取為30，集合卡爾曼濾波的集合數(shù)為30，系統(tǒng)狀態(tài)初始值(x0,y0,z0)仍取為(1.508870，-1.531271，25.46091)，試驗(yàn)中保持模型積分步長(zhǎng)h=0.025，積分區(qū)間為[0,10]不變，依次抽取0，5，10，…，400h時(shí)刻的狀態(tài)量作為N=80個(gè)觀測(cè)數(shù)據(jù)，試驗(yàn)中仍以均值為0、方差為1的高斯隨機(jī)數(shù)作為觀測(cè)噪聲，圖6給出了每隔5步同化一次觀測(cè)的試驗(yàn)中前400步的結(jié)果.

圖6　兩種濾波同化效果比較

從圖6可以看出，兩種濾波方法的狀態(tài)變量曲線都接近于真值曲線，混沌現(xiàn)象基本相同，表明兩種濾波方法都能很好地校準(zhǔn)模型的運(yùn)行軌跡.此外，從圖6還可以看出，在Lorenz混沌系統(tǒng)圖左側(cè)部分，真值曲線與兩種濾波的狀態(tài)變量曲線基本處于重合狀態(tài)；在Lorenz混沌系統(tǒng)圖右側(cè)部分，盡管兩條曲線都十分接近真值曲線，但仍可清晰看出，降秩平方根濾波的狀態(tài)變量曲線處于集合卡爾曼濾波狀態(tài)變量曲線與真值曲線之間，即降秩平方根濾波的狀態(tài)變量估計(jì)值更接近于真值.因此，降秩平方根濾波方法的同化性能在一定程度上可能優(yōu)于集合卡爾曼濾波的同化性能.為進(jìn)一步說明濾波性能，實(shí)驗(yàn)比較了不同濾波情況下，混沌系統(tǒng)三個(gè)狀態(tài)變量x，y，z的均方根誤差.圖7給出了兩種濾波情況下，混沌系統(tǒng)3個(gè)狀態(tài)變量的一維同化效果比較和RMSE值變化.

圖7　兩種濾波同化效果對(duì)比

從圖中可以看出，兩種濾波方法整體較為平穩(wěn)，局部均存在異常波動(dòng).降秩平方根濾波的波動(dòng)幅度明顯小于集合卡爾曼濾波，即短暫發(fā)散程度小于集合卡爾曼濾波，降秩平方根濾波的RMSE值小于集合卡爾曼濾波.由此可以看出，降秩平方根濾波用于三變量Lorenz混沌系統(tǒng)狀態(tài)預(yù)報(bào)時(shí)，同化結(jié)果要優(yōu)于集合卡爾曼濾波.

3結(jié)論

基于Cholesky分解方法提出了降秩平方根濾波算法，并在Lorenz混沌系統(tǒng)上檢驗(yàn)了降秩平方根濾波的性能.同化試驗(yàn)結(jié)果表明，經(jīng)Cholesky方法分解后的卡爾曼濾波在非線性系統(tǒng)中可以避免對(duì)模型算子和觀測(cè)算子進(jìn)行線性化，能夠準(zhǔn)確地更新模式預(yù)報(bào).通過敏感性試驗(yàn)，研究了觀測(cè)方差、同化窗口長(zhǎng)度對(duì)同化效果的影響，比較了降秩平方根濾波和集合卡爾曼濾波的同化效果.

1)試驗(yàn)結(jié)果表明，經(jīng)Cholesky分解的卡爾曼濾波融合模擬觀測(cè)信息能夠?qū)orenz混沌系統(tǒng)狀態(tài)變量進(jìn)行實(shí)時(shí)更新，分解后的卡爾曼濾波同化性能穩(wěn)定.

2)觀測(cè)方差越小，系統(tǒng)同化效果越好，觀測(cè)方差過大容易導(dǎo)致降秩平方根濾波發(fā)散.

3)文中討論了同化窗口長(zhǎng)度對(duì)同化效果的影響.結(jié)果表明，模型運(yùn)行400步以內(nèi)，隨著同化窗口長(zhǎng)度的增大，狀態(tài)變量的均方根誤差值會(huì)隨之增大，說明在同化系統(tǒng)中及時(shí)、足量的觀測(cè)信息對(duì)改善同化效果有很大的作用.

4)與集合卡爾曼濾波相比，分解后的降秩平方根濾波收斂速度明顯提高，在系統(tǒng)短期運(yùn)行時(shí)間內(nèi)(文中為200步)，降秩平方根濾波的同化性能優(yōu)于集合卡爾曼濾波，并且比集合卡爾曼濾波更穩(wěn)定.

文中基于模型運(yùn)行的前400步試驗(yàn)結(jié)果對(duì)兩種同化方法的優(yōu)劣進(jìn)行了分析、比較.試驗(yàn)結(jié)果表明,混沌系統(tǒng)的非線性現(xiàn)象對(duì)同化方法有較大影響，雖與傳統(tǒng)混沌系統(tǒng)研究的結(jié)論相一致，但與傳統(tǒng)混沌系統(tǒng)研究相比，文中模型運(yùn)行時(shí)間較短，模型長(zhǎng)時(shí)間運(yùn)行的情況需進(jìn)一步討論.此外，文中試驗(yàn)及上述結(jié)論均是根據(jù)簡(jiǎn)單數(shù)值試驗(yàn)結(jié)果得到，而降秩平方根濾波應(yīng)用到大氣、海洋、陸面模型中還需作進(jìn)一步研究和討論.

參考文獻(xiàn):

[1]LORENZEN.Deterministicnonperiodicflow[J].Journal of the Atmospheric Sciences，1963，20：130.

[2]劉福才，梁曉明.Hénon混沌系統(tǒng)的廣義預(yù)測(cè)控制與同步快速算法[J].物理學(xué)報(bào)，2005，54(10)：4584.

[3]李麗香，彭海朋，盧輝斌，等.Hénon混沌系統(tǒng)的追蹤控制與同步[J].物理學(xué)報(bào)，2001，50(4)：629.

[4]高山紅，吳增茂，謝紅琴.Kalman濾波在氣象數(shù)據(jù)同化中的發(fā)展與應(yīng)用[J].地球科學(xué)進(jìn)展，2000，15(5)：571.

[5]KALMANRE.Anewapproachtolinearfilteringandpredictionproblems[J].Transaction of the ASME-Journal of Basic Engineering SeriesD，1960，82(2)：35.

[6]EVENSENG.Sequentialdataassimilationwithanonlinearquasi-geostrophicmodelusingMonteCarlomethodstoforecasterrorstatistics[J].Journal of Geophysical Research，1994，99(5)：10143.

[7]SAKOVP，OKEPR.Adeterministicformulationoftheensemblekalmanfilters:analternativetoensemblesquarerootfilters[J].Tellus SeriesA，2008，60(2)：361.

[8]JEFFREYS，WHITAKER，THOMASM，etal.Ensembledataassimilationwithoutperturbedobservations[J].Monthly Weather Review，2002，130(7)：1913.

[9]李昊睿，張述文，邱崇踐，等.四維變分同化和集合平方根濾波聯(lián)合反演土壤濕度廓線的研究[J].大氣科學(xué)，2010，34(1)：193.

[10]王世璋，閔錦忠，陳杰，等.迭代集合平方根濾波在風(fēng)暴尺度資料同化中的應(yīng)用[J].大氣科學(xué)，2013，37(3)：563.

[11]曹小群，宋君強(qiáng)，張衛(wèi)民，等.基于MCMC方法的Lorenz混沌系統(tǒng)的參數(shù)估計(jì)[J].國(guó)防科技大學(xué)學(xué)報(bào)，2010，32(2)：68.

[12]薛具奎.非線性物理基礎(chǔ)[M].蘭州：蘭州大學(xué)出版社，2002.

(責(zé)任編輯孫對(duì)兄)