1、第第7 7講講 大氣資料的四維變分同化方法大氣資料的四維變分同化方法 1 1 四維變分同化基本原理四維變分同化基本原理 四維同化的概念-利用模式消化吸收多時刻觀測，不斷改進預報，優化大氣狀態的估計。 變分方法是一種實現四維同化的有力工具。 在進行大氣資料分析時，我們有兩種基本的可用信息：（1）觀測；（2）大氣遵循的物理規律。前面我們在作資料分析中用到過一些簡化的物理約束，四維變分同化利用完整的大氣模式來作為物理約束。 四維變分同化的基本思想是調整初始場，使由此產生的預報在一定時間區間（同化窗口）內與觀測場距離最小 4DVAR示意圖：按照這樣的思想，四維同化變分同化可以表述為極小化下面的目標泛函

2、：這里x0=x(0), xt=x(t). xt是由下面的預報模式產生的解:離散形式（n=0 成為三維同化）OB101JJ )() )(2121)(dtxHyxHJtttTtBOTBOOOyxxBxxxt y)(xO)y)(x(21xxBxx21xr0r1Tr0rn0rB01TB00H)J(H0)(0),(xoxtxFtX 4DVAR4DVAR是微分方程反問題是微分方程反問題 將已知微分方程和定解條件（初條件，邊條件）求方程的解的問題作為正問題，那末，已知方程的解（部分解）或解的某種函數反求定解條件或者方程的一些未知項的問題被稱之為微分方程的反問題。因此，四維變分同化也是一類微分方程的反問題。

3、求反問題的解的過程稱為反演。我們可將觀測y y近似看作預報模式（方程）的解的某種函數，那末上面表述的四維變分同化就是由觀測反演初值的問題。四維變分同化的一個顯著特點是利用了過去時間的觀測資料，而且同化后的場是模式的一個預報場，不會出現不協調的問題。四維變分同化方法還有能力從一部分觀測變量去反演另外的變量。比如，由高度的觀測反演風場。 關于反問題的進一步討論關于反問題的進一步討論 如果將由“原因”推得“結果”的問題稱為正問題，則由“結果”推求“原因”的問題可稱為反問題。 正問題 z=R(u) 這里算子R R 已知，由u求z為正問題。反問題是原來的已知條件未知（或部分未知），而原問題的解已知，即由

4、z求u的問題，形式上可寫為 u=R-1z 這里R R-1是R R 的逆算子。我們要研究的反問題一般是指R R-1的顯式表達式不可知的情況，只能由u的“表現”z間接推求u。 一個簡單的例子一個簡單的例子 擴散輸送問題擴散輸送問題定解條件：定解條件： 幾類反問題：幾類反問題： 待定微分方程中的未知參數的反問題待定微分方程中的未知參數的反問題算子識別；算子識別； 待定初始條件的反問題待定初始條件的反問題逆時間過程問題；逆時間過程問題； 待定邊界條件的反問題待定邊界條件的反問題邊界控制問題；邊界控制問題； 待定邊界形狀的反問題待定邊界形狀的反問題幾何反問題。幾何反問題。 還有的反問題是幾類相混合。還有

5、的反問題是幾類相混合。0)()(xcxkxxcxutc,baxxx)(), 0(0 xcxc)(),(1tgxtca)(),(2tgxtcb解反問題的主要困難解反問題的主要困難不適定性不適定性 解的存在性、唯一性和穩定性不滿足解的存在性、唯一性和穩定性不滿足 吉洪諾夫的論著不適定問題的解法首先引入“條件適定條件適定”的概念，基本思想是：放棄求精確解轉而求近似解解決了解不存在的困難；近似解總存在，但不唯一，此時再加適當約束條件，找出具有穩定性的解來。 廣義解，目標泛函廣義解，目標泛函 反問題Au =z 的廣義解： uU， zZ 對于給定的zZ在集U上使 取極小值的u*U, 稱為方程Au=z在U上

6、的廣義解。 若zAU，廣義解等于經典解。(AU為U的映象), 為距離。經典解解不存在: z不屬于AUzAu . 假定算子方程 Au =z 的逆算子A-1存在但不連續依賴于z，由u=A-1z計算u不再現實。正則化的思想是構造一個連續的算子去逼近A-1，從而得到穩定的（但是近似的）解。具體而言，他將求穩定的反問題的解歸結為求下面泛函的極小值： 非負泛函,:觀測誤差.正則系數，比如罰函數：)(,uzAuuzM)(udxdxdYHWIbaxx2222)(如何定? 給出讓 迭代。zAu 2, 2, 最優控制理論最優控制理論 四維變分同化將問題提為一個最優控制問題。四維變分同化將問題提為一個最優控制問題。

7、最優控制問題的一般提法：最優控制問題的一般提法： 在系統工程中，一個系統通常可以用在系統工程中，一個系統通常可以用n n個變量完個變量完全描述清楚。我們以向量全描述清楚。我們以向量X X 記這記這n n個變量，并稱之個變量，并稱之為狀態向量。系統的運動方程可以用時間間隔為狀態向量。系統的運動方程可以用時間間隔T T1 1,T,T2 2上的狀態方程來表示：上的狀態方程來表示：,),(21TTttuxfX 其中其中是是n 維狀態向量；維狀態向量；u是是r維（維（r n）控制向量，）控制向量，它是從系統之外按一定要求施加到系統上來的；它是從系統之外按一定要求施加到系統上來的；f是是n 維向量函數。維

8、向量函數。由方程可以看出，只要由方程可以看出，只要f 滿足一定的條件，在確定滿足一定的條件，在確定的初始狀態下，如果在時間間隔的初始狀態下，如果在時間間隔t0, TT1, T2上給定了一個控制規律上給定了一個控制規律uu（t）, 那末狀態方程在那末狀態方程在t0,T上將有唯一解。這個解表示了系統的相點上將有唯一解。這個解表示了系統的相點在在n維狀態空間的一個運動維狀態空間的一個運動 控制規律控制規律u不同時，相應的系統運動也是不不同時，相應的系統運動也是不同的所謂最優控制問題就是要選擇適當的同的所謂最優控制問題就是要選擇適當的控制控制 規律規律 u(t) 使相應的品質指標使相應的品質指標(co

9、st function, 代代價函數或目標函數價函數或目標函數) 達到最小這里第一項代表了對終端條件的達到最小這里第一項代表了對終端條件的要求，第二項代表了對整個變化過程的要求，要求，第二項代表了對整個變化過程的要求，而總的品質指標是這兩方面要求的綜合而總的品質指標是這兩方面要求的綜合 rtdttuxLTTxKuJ0),(),()( 最優控制問題的一般提法應為：給定了系統的狀態最優控制問題的一般提法應為：給定了系統的狀態方程以及品質指標和初始條件（方程以及品質指標和初始條件（t0, x0），目標集），目標集，控制域控制域后，要尋找一個容許控制后，要尋找一個容許控制u*（t），使系統），使系統在

10、時間間隔在時間間隔t0, T上由初始狀態上由初始狀態X0轉移到目標集轉移到目標集上的某一點上的某一點Xr，且使相應的品質指標，且使相應的品質指標J（u）為極）為極小這里終端時間小這里終端時間T可以是固定的，也可以是自由可以是固定的，也可以是自由的的 (容許控制的作用下，系統由初始狀態容許控制的作用下，系統由初始狀態X0轉移到終端轉移到終端狀態狀態X有些情況，我們要求終端狀態應該是狀態有些情況，我們要求終端狀態應該是狀態空間中某一個點集空間中某一個點集中的一個點集合中的一個點集合常稱為目常稱為目標集當標集當只是狀態空間中的一個確定的點時，稱只是狀態空間中的一個確定的點時，稱終端是固定的當終端是固

11、定的當是整個狀態空間時，稱終端是是整個狀態空間時，稱終端是自由的自由的) 根據這些要求所求得的控制函數稱為最優控制或極值控制，記為u*（t）而相應的狀態向量稱為最優軌跡或極值值軌跡，記為X*(t） 從本質上講，所求得的最優控制從本質上講，所求得的最優控制u*是時間是時間的函數但從形式上來看，它可以直接表示的函數但從形式上來看，它可以直接表示為時間的函數為時間的函數u*（t），也可以表示為狀態），也可以表示為狀態的函數的函數u*（x）如果求得的最優控制被直）如果求得的最優控制被直接表示為時間的函數，即接表示為時間的函數，即u*u*（t），則），則稱所求得的是開環最優控制若求出的最優稱所求得的是開

12、環最優控制若求出的最優控制被表示為狀態的函數，即控制被表示為狀態的函數，即 u*u*(x）則稱所求得的是閉環最優控制從控制理論則稱所求得的是閉環最優控制從控制理論的一般知識可以知道，閉環控制比開環控制的一般知識可以知道，閉環控制比開環控制有不少優越之處，所以在實際應用中，總是有不少優越之處，所以在實際應用中，總是力圖尋找閉環控制規律但是，在實際計算力圖尋找閉環控制規律但是，在實際計算中，尋找閉環最優控制規律要比尋找開環最中，尋找閉環最優控制規律要比尋找開環最優控制規律困難得多優控制規律困難得多 這里所謂“最優”的標準就是品質指標J 取極小。但J 的形式的確定完全取決于對本身的具體要求以及工作人

13、員的工作經驗，并無普遍適用的法則。可以說，許多最優控制問題的最困難之點就在于無法提出一個恰當的品質指標泛函。 3, 3, 用變分法解最優控制問題用變分法解最優控制問題 顯然最優控制問題是一個求帶約束的泛函顯然最優控制問題是一個求帶約束的泛函的極小問題，前面講述的變分方法同樣適的極小問題，前面講述的變分方法同樣適用于解決這一問題。以下用一個例子來說用于解決這一問題。以下用一個例子來說明一般做法。明一般做法。0)()(xcxkxxcxutc反演初值和反演初值和k(x), 定義目標泛函定義目標泛函：babaxxBxxobsTdxxkxkdxdtccxtwJ)()()(,(21220baxxTdxdt

14、xcxkxxcxutcxtcJcckL)()(),(),(*0*0無約束極小：無約束極小：dxdtxcxkxxcxutccdxdtxcxkxxcxkxxcxutcxtcJLbabaxxTxxT)()()()()(),(*0*0變分：變分：分部積分得到：babaxxBxxkobscTkdx)kk(wcdxdt)cc(wL0dxdt)xc)x(k(xxc)x(utccdxdt)xccxk)xc)x(k(xxc )x(utcc(babaxx*Txx*T00baxxdx)x ,(*c )x ,( c00利用了邊界： 0)x , t ( c)x , t ( cba0),(*),(* , 0),(*bxt

15、cxtcxTca要求終端和邊界：要求變分為零得到伴隨方程要求變分為零得到伴隨方程（它的定解條件前面已經（它的定解條件前面已經給出）給出）兩點邊值問題兩點邊值問題 ： （迭代求解）迭代求解）)()()(*obscccwxcxkxxcxutc0)()(xcxkxxcxutc0*)(0TBkkdtxccxkkwJ用下降算法用下降算法, ,梯度：梯度：)x ,(*cJc00TBkkdtxc*cx)kk(wJ0 這里這里x x0 0= =x x(0), (0), x xt t= =x x(t)(t)，x xt t是由下面的預報模式產生的是由下面的預報模式產生的解解: :控制變量也可以是方程的參數，比如控

16、制變量也可以是方程的參數，比如s(x)s(x)OB101JJ )() )(2121)(dtHOHBJttTTOxyxyxxxxxtttBOBO0)(0),(xxxxotFt回到四維變分同化來，也是一個最優控制問題極回到四維變分同化來，也是一個最優控制問題極小化下面的目標泛函：小化下面的目標泛函：)()(xxxsFt 從另一角度來推導：先對伴隨算子做些說明。在希爾伯特（Hilbert）空間（H空間），稱A*為線性算子線性算子A A 的伴隨算子，只要對于屬于H的任意元素x，y均有 =, 為內積. 在Hilbert空間，內積的定義： (在N 維向量空間內積的定義顯然，N 維向量空間的矩陣算子A的伴隨

17、算子就是其轉置AT ，因為xTAy= yT ATx)xydyx,NIiiTyxyxyx1,預報模式方程：線性化得到切線性模式: 伴隨模式:0F(u)ut0)(0,uuuFuoTttu（7.2.1)(7.2.2)0),(),(*TxttuuguFu（7.2.3)F F*是 的伴隨算子(轉置）：將u u*和(7.2.2)作內積，減去uu和(7.2.3)作內積,然后在(0,)區間對時間積分得到（空間邊界變分為零條件下）:uFt 算子 的伴隨算子t dttTTTT),(00*ugux)(u* *uFuuFuTuT如果目標函數是：讓得到表明J對初值的梯度是 。dtLdtJTTO),(21(21)(00t

18、ttt1tTttuy)H(uyO)H(uyu tt1tTuHyOuHuugtLt),( JdtLdt000uuuHyOxHuu)(Ttt1tTTT)0(4 4 四維變分同化的計算過程四維變分同化的計算過程實際問題中，微分方程都被離散化，這時目標函數實際問題中，微分方程都被離散化，這時目標函數為為控制方程：控制方程：這里這里n n是有觀測的總時間層數，將目標函數是有觀測的總時間層數，將目標函數) )寫為寫為 )()(2121)(100100rrrTrrrnrBTByxHOyxHxxBxxxJ0)()(xxQtxxrrrronrbobJJJJJ)(0（7.4.1)為了獲得x0的最優估計x*0, 需

19、計算J 相對x0的梯度：其中 的計算需要對伴隨模式積分。由于x0的擾動x0引起的J0的擾動由（7.4.1）略去高階小量,由上面兩式得到obJJJ)(1BbxxBJoJ20000000)()()(xOxxJxJxxJxJToooo2100)()(rrTrrrrTrnroxOxyxHOHxJrTrrrrTrnrToxyxHOHxxJ1000（9.4.6）對預報模式作擾動得到切線性模式：代入(7.4.6)中的 , (7.4.6)取極限 ，(7.4.6)意味著 (7.4.7)記為 0)()(xxQtxxrrr 0 xxPtxrrrx rTrrrrTrnrToxyxHOHxxJ1000rrrrTrTrR

20、royxHOHPxJ)()(10000 xrRroxxJ000)(即下面方程的解：rx0)( )(, 1 , 0,)()( 01rTrrrrrrTrrtxxPxRryxHOHtx(7.4.9)(7.4.8)即每一個 是由時間 出發，以 為初值向后積分伴隨方程到t=t0得到。由于方程（7.4.8-9）都是線性的，實際上只須對伴隨模式由tr到t0作一次積分，初值為零，而在到達觀測時刻tr時，伴隨變量加上如(7.4.9)右端表示的一個強迫項。rx0rtrrrrTryxHOH)(1RryxHOHtxtxxPxxxPPPxPPxPxxxPPPxPPxxPxrrrrTrrrTrrRRTRTTTTTRRRR

21、RRRRR, 1 , 0,)()( )( )( P.P .)(10RT1102101000R0021221111實際過程實際過程舉例：預報X(1,0)=a，X(2,0)=b，X(1,1)= X(1,0)* X(2,0)+c X(1,0)X(2,1)= X(1,0) +d X(2,0)X(1,2)= X(1,1)* X(2,1)+c X(1,0)X(2,2)= X(1,1) +d X(2,1)目標函數：J=0.5（X(1,2)-y1）*2+( X(2,2) y2)*2線性方程X(1,0)= a，X(2,0)= b， X(1,1)= X(2,0)* X(1,0)+ X(1,0)* X(2,0)+c

22、 X(1,0) X(2,1)= X(1,0)+ dX(2,0）0020101020211111xpxxxdxcxxx0211121111121221221xppxpxxxdxcxxx需要求222112221222122/*yxyxxJxJxxx)/)(/()/)(/()/)(/()/)(/(/*211212212222112222111212211121111xxxJxxxJxxxJxxxJxJxJxxx2122121*/xPPTTxJxJ210100*xPPxPxTTT201020100/*xJxJxxx先計算然后類似：伴隨碼（伴隨碼（Adjoint code)生成技術生成技術 為了得到目標

23、函數對初值的梯度我們必須積為了得到目標函數對初值的梯度我們必須積 分伴隨模式。伴隨模式方程實際是原模式的切線分伴隨模式。伴隨模式方程實際是原模式的切線 性方程的伴隨方程。伴隨模式的解析形式只能作性方程的伴隨方程。伴隨模式的解析形式只能作 為理論推導用，實際問題是離散化的，預報程序為理論推導用，實際問題是離散化的，預報程序 中還有些是不能寫成解析公式的。要保證相應的中還有些是不能寫成解析公式的。要保證相應的 伴隨模式嚴格成立，通常的作法是先根據原模式伴隨模式嚴格成立，通常的作法是先根據原模式 計算程序寫出切線性模式程序，再直接根據切線計算程序寫出切線性模式程序，再直接根據切線 性模式程序一一對應

24、地寫出伴隨程序。一個預報性模式程序一一對應地寫出伴隨程序。一個預報 模式的程序有上萬條語句，首先寫出他的切線性模式的程序有上萬條語句，首先寫出他的切線性 模式程序，然后根據切線性模式程序寫出伴隨程模式程序，然后根據切線性模式程序寫出伴隨程 序，工作量是巨大的。序，工作量是巨大的。這里我們用一些簡單例子介紹一般規則。這里我們用一些簡單例子介紹一般規則。在預報模式中，一個循環過程一般可以表示在預報模式中，一個循環過程一般可以表示為矩陣與向量的乘為矩陣與向量的乘Y YAXAX這里這里A A與與X X無關（線性），則切線性程序和原來相同，無關（線性），則切線性程序和原來相同，Y YAXAX伴隨過程為伴

25、隨過程為X X* *A AT TY Y* *如果如果A A與與X X有關（非線性），那末切線性程序為有關（非線性），那末切線性程序為Y YA AX X伴隨過程為伴隨過程為X X* *A AT TY Y* *例如 DO I=1,N-1 X(I)= a * Y(I+1)*2 END DO 切線性程序為: DO I=1,N-1 X(I)= a * 2*Y0(I+1)* Y (I+1) END DO為了以下方便，記a * 2*Y0(I+1)B（I）顯然，Y是輸入變量，X是輸出變量。注意，這里Y0是“基本變量”，X和Y是擾動變量。基本變量是由原模式產生的。 需要區分以下兩種情況：需要區分以下兩種情況：（1 1）只有）只有X X是輸出變量，是輸出變量，Y Y在以后不在用到（在以后不在用到（Y Y不是不是輸出變量）；輸出變量）；（2 2）X X和和Y Y都是輸出變量，就是說都是輸出變量，就是說Y Y在以后還要到。在以后還要到。對第一種情況，上述循環的矩陣是對第一種情況，上述循環的矩陣是 )(00000)3(000000)2(01NBBBA(N -1)N對第二種情況，矩陣是 (2N -1)NNBBBA000000000000000010000011000010000