同化算法課程

同化1.基本概念1.1分析在氣象學(xué)及地球物理的其他分支中確定某時(shí)刻的物理系統(tǒng)得真實(shí)狀態(tài)的過(guò)程稱(chēng)為分析。作為分析的基礎(chǔ)的信息包括觀測(cè)數(shù)據(jù)以及物理系統(tǒng)及其初邊值條件，還可能有對(duì)分析的約束條件。分析本身對(duì)描述物理系統(tǒng)是有用的，它也可以作為研究系統(tǒng)演化的初始狀態(tài)。例1：經(jīng)典分析方法。在計(jì)算機(jī)的發(fā)明之前，氣象學(xué)家是用手工做氣象分析的。觀測(cè)是在大致同時(shí)（例如在船上報(bào)告12UTC（CoordinatedUniversalTime協(xié)調(diào)世界時(shí)）表面壓力）對(duì)溫度，風(fēng)及壓力作記錄。這些手工的分析比僅僅用觀測(cè)數(shù)據(jù)作空間的插值得到的更多，有些人把它用于早期數(shù)值天氣預(yù)報(bào)。方法如下：把預(yù)報(bào)模型在格點(diǎn)處的首次猜測(cè)場(chǎng)插值到觀測(cè)點(diǎn)，把觀測(cè)值與插值的差插值回格點(diǎn)作為訂正（隨著格點(diǎn)與觀測(cè)點(diǎn)的距離的增大而減少）。這個(gè)方法的重要的一步就是在對(duì)觀測(cè)值的插值時(shí)首次猜測(cè)場(chǎng)的引入。雖然這個(gè)早期的工作引出了最優(yōu)插值法，但它仍然是廣泛使用的分析方法。例2：臭氧分析。臭氧是在空間用遙感測(cè)量的。一天得到的數(shù)據(jù)覆蓋了大部分地球，如果用分析的方法來(lái)把空缺部分的數(shù)據(jù)補(bǔ)充完整從而得到全球的臭氧圖。但是數(shù)據(jù)是在24小時(shí)的時(shí)間區(qū)間上得到的，而臭氧場(chǎng)在某些地方某些時(shí)刻可能有很大的變化，所以我們?cè)诳臻g的插值效果可能很差。一個(gè)更好的分析方法是把臭氧的動(dòng)態(tài)演化和觀測(cè)的時(shí)間考慮在內(nèi)。這樣就可以得到符合大氣臭氧演變的全景圖。例3．天氣預(yù)報(bào)：分析者的前期準(zhǔn)備工作是：以初始條件為天氣預(yù)報(bào)的基礎(chǔ)，加上新得到的觀測(cè)值、把以前分析得到的預(yù)報(bào)值作為大氣狀態(tài)的背景場(chǎng)的估計(jì)值。在分析的準(zhǔn)備工作中，一個(gè)特別要做的工作通常是用附加的約束條件來(lái)衰減重力波。在地球物理系統(tǒng)的時(shí)間演化的預(yù)報(bào)中一個(gè)特別的挑戰(zhàn)就是系統(tǒng)的非線性性及對(duì)初始條件的敏感性。我們熟知，即使對(duì)于一個(gè)完美的預(yù)報(bào)模型，任何預(yù)報(bào)在有限時(shí)間后就會(huì)和真實(shí)狀態(tài)發(fā)生很大的偏離。這個(gè)例子就是Duffing方程：假設(shè)真實(shí)狀態(tài)是X+0.05X+x3二7.5cost,x（0）二3,x（0）二4而把分析誤差解釋為初始條件，分析的初始條件為x（0）二3.01,X（0）二4.01，則在時(shí)間t二35以后解就發(fā)生越來(lái)越大的偏離。1.2同化如上討論的，一個(gè)分析可以很簡(jiǎn)單，比如空間的插值。然而，更好的結(jié)果可以通過(guò)在分析中考慮到物理系統(tǒng)的動(dòng)態(tài)演化。一個(gè)結(jié)合了關(guān)于時(shí)間分布的觀測(cè)值以及動(dòng)態(tài)模型的分析方法稱(chēng)為同化或數(shù)據(jù)同化。同化問(wèn)題可以從多角度來(lái)討論，取決于背景及喜好（比如控制論，估計(jì)理論，概率論，變分分析等等）。在數(shù)值模式中的算法中最常用的是離散方法。然而，在此介紹中要說(shuō)明的是不同的同化算法是怎樣從同一個(gè)常用來(lái)源得到的。并對(duì)每種算法給出其近似程度（特別是最優(yōu)插值，三維變分，四維變分，卡爾曼濾波）。其中最吸引人的方法是從變分的觀點(diǎn)討論連續(xù)的同化問(wèn)題。先來(lái)看同化問(wèn)題和熟悉的初值問(wèn)題的不同之處。同化問(wèn)題可以看成帶有誤差的初值問(wèn)題。模式狀態(tài)x的方程帶有模式誤差F，初始值帶有背景誤差F，觀測(cè)值y帶有觀mbn測(cè)誤差￡。o,n+M(x)="e-<Qtm,y=H(x)+s,n=nn0,nx(0)=x(0)+rbb其中H是觀測(cè)算子，通過(guò)模式計(jì)算得出格點(diǎn)n處的觀測(cè)值y。

nn1.3以概率方式的描述1.3.1概率密度函數(shù)上述同化問(wèn)題含有未知的誤差，故顯然其解是統(tǒng)計(jì)意義上的。因此需要了解有關(guān)誤差的統(tǒng)計(jì)特點(diǎn)。對(duì)于特定系統(tǒng)的觀測(cè)誤差和模式誤差的特征的刻畫(huà)是數(shù)據(jù)同化的另一個(gè)基礎(chǔ)。最理想的情況是得到誤差的全概率密度函數(shù)P($,巧,巧,t)，但實(shí)際上我們對(duì)P了解甚少，至多知mbo道其期望值及協(xié)方差。在大部分情況下，我們假定模式誤差、背景誤差和觀測(cè)誤差是互不相關(guān)的。這時(shí)，全概率密度函數(shù)是各密度函數(shù)(模式誤差概率密度函數(shù)P，背景誤差概率密度函數(shù)P，觀測(cè)誤mb差概率密度函數(shù)P)的乘積：P二PPP二exp(lnP+lnP+lnP)。ombombo看起來(lái)最佳的分析解對(duì)應(yīng)于P的最大值。但實(shí)際上不總是這樣的。例：考慮一個(gè)一維的概率P,設(shè)在觀測(cè)值為y時(shí)發(fā)生x的條件概率為P(x1y)，對(duì)于一個(gè)高斯分布，x的最大似然點(diǎn)和x的均值點(diǎn)是重合的，故選取密度函數(shù)的最大值是最佳的。但對(duì)于一個(gè)非高斯分布，x的最大似然點(diǎn)是不同于x的均值點(diǎn)。這時(shí)x的最大似然點(diǎn)并不是一個(gè)最好的解。因此我們必須要研究實(shí)際的密度函數(shù)以保證最大似然解的確是我們要的分析解。1.3.2代價(jià)函數(shù)概率的最大化常常表示為概率密度的負(fù)指數(shù)的最小化：(稱(chēng)之為代價(jià)函數(shù)J)minJ二min(-lnP—lnP—lnP)二min(J+J+J)mb0mbo例:假設(shè)忽略模式誤差，且P=1exp(—FB-￡/2)，P=1exp(—FR-€/2)，b耳；2兀IBlbbo存|R|oo其中B是背景誤差協(xié)方差陣，R是觀測(cè)誤差協(xié)方差陣，則忽略不影響結(jié)果的常數(shù)后代價(jià)函數(shù)為J=(FB-1蘆+FR-1蘆)/2bboo=((X(0)—x(0))'B-1((x(0)—x(0))+(y—H(X))'R-i(y—H(x)))/2bb當(dāng)觀測(cè)算子是線性算子時(shí)，代價(jià)函數(shù)是x的二次函數(shù)，這時(shí)有惟一的最小解。而且最小解的計(jì)算是高效的。但在實(shí)際問(wèn)題中，模式和觀測(cè)算子都是非線性的。得到的只能保證是局部最小解，不一定是全局最小解。不過(guò)，對(duì)于高斯誤差，有很有效的算法，它對(duì)代價(jià)函數(shù)作局部的線性化疊代而得到最優(yōu)解。在氣象學(xué)中，高斯分布通常是一個(gè)很好的近似分布。2．變分?jǐn)?shù)據(jù)同化誤差協(xié)方差陣與偏差協(xié)方差陣

模式誤差協(xié)方差陣Q=Q(g,t,了,t')描述在兩個(gè)不同相時(shí)點(diǎn)的協(xié)方差，階數(shù)同方程中未知量的個(gè)數(shù)。背景誤差協(xié)方差陣B=B(C,◎)描述在兩個(gè)不同相點(diǎn)的協(xié)方差觀測(cè)誤差協(xié)方差陣R分析誤差協(xié)方差陣A=A(C,了)，2.1.2偏差偏差是同化中的一大困難，在協(xié)方差陣中可能會(huì)有非零的偏差，而實(shí)際上卻很難確定它。比如不能確定偏差是由觀測(cè)還是模型預(yù)報(bào)來(lái)的。在同化算法中怎樣加入偏差以及如何估計(jì)偏差都是具有挑戰(zhàn)性的工作。一旦知道的偏差，我們可以通過(guò)減去這個(gè)偏差而得到一個(gè)無(wú)偏的量，這稱(chēng)為無(wú)偏化。2.2變分同化問(wèn)題的描述對(duì)于誤差是高斯分布的情況，最大似然解就是最佳的解，這時(shí)代價(jià)函數(shù)是二次型：J(x)二1任,t)Q-G(g,tr)d?d&dtdt'2mm+1任)B-G(C)d?dC+1k'(g,t)R-代(了,t，)dgd^'dtdt'2bb2oo其中關(guān)于時(shí)間的積分是在同化窗口［0,］上進(jìn)行的，f(C,t)={^(C,t)(y-H(X))},onnnTOC\o"1-5"\h\zP(C,t)是一局部化函數(shù)：W(C,t)dCdt二1nn設(shè)X是J(X)的最小解，令x(C,t)=X(C,t)+Yn(C,t)，則對(duì)于任意數(shù)丫及任意函數(shù)aan(C,t)，g(Y)：=J(X(C,t))一J(x(C,t))、0，故—0，即adtY-0dg(Y)dtpY-0、k(C,0)B-1(X(C',0)眄dMJ門(mén)dtdx丿a-Xb(C,0))dCdCbdX(C,t')+M(X(C',t'))dC'dt'丿adta-J［型吒］R-i(y—H(X))dCdC'dtdt'=0ldx丿a記兀(C,t):—JQ-1■adt血(C,t')+M(記兀(C,t):—JQ-1■adt分部積分得Jn(g冗(g,t)dg+