擺玉龍，尤元紅，邵　宇，張轉花

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

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

擺玉龍，尤元紅，邵宇，張轉花

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

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

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

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

E-mail:yulongbai@gmail.com

基金項目:國家自然科學基金資助項目(41461078)

作者簡介：擺玉龍(1973—)，男，甘肅會寧人，教授，碩士研究生導師.主要研究方向為數(shù)據(jù)同化和參數(shù)估計.

中圖分類號：TP 309.7

文獻標志碼：標志碼：A

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

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

1理論方法

1.1卡爾曼濾波

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

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

2)計算增益矩陣.利用i時刻的預報誤差協(xié)方差矩陣計算增益矩陣.

其中H為觀測矩陣；HT為觀測矩陣的轉置.

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

1.2降秩平方根濾波

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

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

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

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

3)分析.

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

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

2數(shù)值試驗

2.1Lorenz混沌系統(tǒng)

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

(16)

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

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

2.2降秩平方根濾波基本同化試驗

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

圖2　第400步同化效果圖

2.3觀測方差對同化效果影響

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

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

圖3　觀測方差對同化效果的影響

2.4觀測窗口長度對同化效果影響

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

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

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

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

圖4　觀測窗口對同化效果影響

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

2.5加密觀測

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

圖5　采用加密觀測前后比較

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

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

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

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

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

圖7　兩種濾波同化效果對比

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

3結論

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

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

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

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

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

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

參考文獻:

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

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

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

[4]高山紅，吳增茂，謝紅琴.Kalman濾波在氣象數(shù)據(jù)同化中的發(fā)展與應用[J].地球科學進展，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].大氣科學，2010，34(1)：193.

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

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

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

(責任編輯孫對兄)