同化算法課程

1、同化1.基本概念1.1分析在氣象學(xué)及地球物理的其他分支中確定某時(shí)刻的物理系統(tǒng)得真實(shí)狀態(tài)的過(guò)程稱為分析。 作為分析的基礎(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（Coordinated Universal Time協(xié)調(diào)世界時(shí)）表面 壓力）對(duì)溫度，風(fēng)及壓力作記錄。這些手工的分析比僅僅用觀測(cè)數(shù)據(jù)作空間的插值得到 的更多，有些人把它用于早期數(shù)值天氣預(yù)報(bào)。方法如下：把預(yù)報(bào)模型在格點(diǎn)處的首次

2、猜 測(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.天氣

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.5cos t, x（0） = 3,x （0） 二 4而把分析誤差解釋為初始條件，分析的初始條件為x（0） = 3.01, X （0

4、） = 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)模型的分 析方法稱為同化或數(shù)據(jù)同化。同化問(wèn)題可以從多角度來(lái)討論，取決于背景及喜好（比如控制論，估計(jì)理論，概率論， 變分分析等等）。在數(shù)值模式中的算法中最常用的是離散方法。然而，在此介紹中要說(shuō) 明的是不同的同化算法是怎樣從同一個(gè)常用來(lái)源得到的。并對(duì)每種算法給出其近似程度 （特別是最優(yōu)插值，三維變分，四維變分，卡爾曼濾波）。其中最吸引人的方法是從變 分的觀點(diǎn)討論連續(xù)的同化問(wèn)題。

5、先來(lái)看同化問(wèn)題和熟悉的初值問(wèn)題的不同之處。同化問(wèn)題可以看成帶有誤差的初值問(wèn)題。模式狀態(tài)x的方程帶有模式誤差廣，初始值帶有背景誤差廣，觀測(cè)值J帶有觀 mbno,n測(cè)誤差8。_初m ，x(0) = x (0) + uVbb其中氣是觀測(cè)算子，1.3以概率方式的描述1.3.1概率密度函數(shù)上述同化問(wèn)題含有未知的誤差，+ M(X) = -_通過(guò)模式計(jì)算得出格點(diǎn)n處的觀測(cè)值yn。y = H (x) + s, n = 1,.,N故顯然其解是統(tǒng)計(jì)意義上的。因此需要了解有關(guān)誤差的統(tǒng)計(jì) 特點(diǎn)。對(duì)于特定系統(tǒng)的觀測(cè)誤差和模式誤差的特征的刻畫(huà)是數(shù)據(jù)同化的另一個(gè)基礎(chǔ)。最理想 的情況是得到誤差的全概率密度函數(shù)P(m，匕,七,

6、t)，但實(shí)際上我們對(duì)？了解甚少，至多知 道其期望值及協(xié)方差。在大部分情況下，我們假定模式誤差、背景誤差和觀測(cè)誤差是互不相關(guān)的。這時(shí)，全概 率密度函數(shù)是各密度函數(shù)(模式誤差概率密度函數(shù)Pm，背景誤差概率密度函數(shù)p，觀測(cè)誤差概率密度函數(shù)P)的乘積：P = P PP = exp(lnP + lnP + lnP )。om b ombo看起來(lái)最佳的分析解對(duì)應(yīng)于P的最大值。但實(shí)際上不總是這樣的。例：考慮一個(gè)一維的概率P，設(shè)在觀測(cè)值為y時(shí)發(fā)生x的條件概率為P(xly)，對(duì)于一個(gè)高 斯分布，x的最大似然點(diǎn)和x的均值點(diǎn)是重合的，故選取密度函數(shù)的最大值是最佳的。但對(duì) 于一個(gè)非高斯分布，x的最大似然點(diǎn)是不同于x的均

7、值點(diǎn)。這時(shí)x的最大似然點(diǎn)并不是一個(gè) 最好的解。因此我們必須要研究實(shí)際的密度函數(shù)以保證最大似然解的確是我們要的分析解。1.3.2 代價(jià)函數(shù)概率的最大化常常表示為概率密度的負(fù)指數(shù)的最小化：(稱之為代價(jià)函數(shù)J)min J = min (- lnP lnP lnP) = min(J + J + J )m b 0m b o例：假設(shè)忽略模式誤差，且P =1exp(-WB-區(qū) /2)，P =1exp(-WR-g /2)，b2兀 | B| b b o：2兀 | R|o o其中B是背景誤差協(xié)方差陣，R是觀測(cè)誤差協(xié)方差陣，則忽略不影響結(jié)果的常數(shù)后代價(jià)函 數(shù)為J = &B-區(qū) +U，R-區(qū))/2=(X(0)-X (

8、0)而-1(眼(0)-x (0) + (y-H保)R-1(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ù)同化2.1誤差協(xié)方差陣與偏差2.1.1協(xié)方差陣模式誤差協(xié)方差陣Q=Q(W t, , t)描述在兩個(gè)不同相時(shí)點(diǎn)的協(xié)方差，階數(shù)同方程中未知量 的個(gè)數(shù)。背景誤差協(xié)方差陣B=B(g,；)描述在兩個(gè)不同

9、相點(diǎn)的協(xié)方差觀測(cè)誤差協(xié)方差陣R分析誤差協(xié)方差陣A=A0頊)，2.1.2偏差偏差是同化中的一大困難，在協(xié)方差陣中可能會(huì)有非零的偏差，而實(shí)際上卻很難確定它。比 如不能確定偏差是由觀測(cè)還是模型預(yù)報(bào)來(lái)的。在同化算法中怎樣加入偏差以及如何估計(jì)偏差 都是具有挑戰(zhàn)性的工作。一旦知道的偏差，我們可以通過(guò)減去這個(gè)偏差而得到一個(gè)無(wú)偏的量，這稱為無(wú)偏化。2.2變分同化問(wèn)題的描述對(duì)于誤差是高斯分布的情況，最大似然解就是最佳的解，這時(shí)代價(jià)函數(shù)是二次型：J(x) = 2芝(；,t)Q-nT (礦tXgdgdtdt+ 2 Jg； (g)B-i% 育)d；d己 + 2廣(如)R- 育,t)dgdM，dtdt 其中關(guān)于時(shí)間的積

10、分是在同化窗口0,T 上進(jìn)行的，E (C,t)=W (C,t)(y H (X)， onn n甲(C,t)是一局部化函數(shù)：m (?,t)d?dt = 1 nn設(shè)X是J(x)的最小解，令X(C,t) = X (C,t) +而(如)，則對(duì)于任意數(shù)Y及任意函數(shù) aan&t)，g(Y)：=J(X(如)-J(Xa(如)0,故譬 =0,即Y=0J dgdt*|聽(tīng)Q-18texJ-i(X (耿0) - x 含,0)d 灑ab-J 竺 1 R-i (y - H(X )dgdgdtdt=0版Jadg(Y)dty=。JiT (；0) B竺 廠所(+M(X &,t，) d?dtJ記兀(?,t) ：=J Q-1成(頊/、版(x (?,t ,)dgdt J分部積分得Jn (& 具)兀(& 具)dgdgdl TOC o 1-5 h z +椅_竺2*血兀(8