條件非線性最優擾動（CNOP）：簡介與數值求解

孫國棟穆穆段晚鎖王強彭飛

（1 中國科學院大氣物理研究所，大氣科學和地球流體力學數值模擬國家重點實驗室，北京 100029；2 復旦大學大氣科學研究院，上海 200433；3 中國科學院海洋研究所，海洋環流與波動重點實驗室，青島 266071；4 中國科學院大學，北京 100049）

條件非線性最優擾動（CNOP）：簡介與數值求解

孫國棟1,4穆穆2段晚鎖1王強3彭飛1,4

（1 中國科學院大氣物理研究所，大氣科學和地球流體力學數值模擬國家重點實驗室，北京 100029；2 復旦大學大氣科學研究院，上海 200433；3 中國科學院海洋研究所，海洋環流與波動重點實驗室，青島 266071；4 中國科學院大學，北京 100049）

介紹了條件非線性最優擾動（Conditional Nonlinear Optimal Perturbation，CNOP）的定義及其在大氣和海洋等可預報性研究中的應用。根據研究對象不同，CNOP分為與初始擾動有關的CNOP（CNOP-I）方法、與模式參數擾動有關的CNOP（CNOP-P）方法和同時考慮初始擾動和模式參數擾動的CNOP方法。目前，CNOP-I方法已經應用于ENSO、黑潮和阻塞可預報性以及熱鹽環流和草原生態系統穩定性的研究。此外，CNOP-I方法也被應用于探討臺風目標觀測的研究，利用CNOP-I方法能夠識別出臺風預報的初值敏感區，通過觀測系統模擬試驗表明在初值敏感區增加觀測能夠有效改進臺風的預報技巧。CNOP-P方法也在ENSO和黑潮可預報性以及熱鹽環流和草原生態系統穩定性研究中得到了應用。為了將CNOP方法應用于更多的領域，本文利用一個簡單的Burgers方程，介紹了如何通過建立Burgers方程的切線性模式和伴隨模式，從而利用非線性最優化算法計算獲得CNOP。這一數值試驗為將CNOP方法應用于更多的領域提供了借鑒。

條件非線性最優擾動方法（CNOP），可預報性，目標觀測

0 引言

大氣和海洋科學中的數值天氣和氣候可預報性研究一直是國內外研究的熱點問題。很多方法已經建立并且應用去探討數值天氣和氣候的可預報性問題。Leith[1]利用蒙特卡羅技術構建大量的隨機初始誤差，研究發現這些初始誤差對6～10d的風場的預報技巧有較大改進。但是隨機樣本個數的選取是非常難以決定的，并且選取的初始誤差是隨機的。近些年，很多學者建立和應用一些方法探討了數值天氣和氣候的可預報性，例如繁殖模向量和集合卡爾曼濾波方法等[2-3]。其中，一些產生動力約束的初始擾動的線性方法被建立去研究數值天氣和氣候可預報性。例如，線性奇異向量方法和Lyapunov向量方法等。這些方法已經被應用于研究海—氣耦合模式中誤差發展的有關規律[4-5]。這些方法屬于初始誤差發展的線性理論[6]。只有在時間較短和初始誤差較小的情況下，線性近似才是成立的。然而，大氣和海洋科學中的物理問題通常是非線性問題。隨著初始誤差的時間演變，將線性誤差理論應用于非線性問題的研究是不合適的。建立非線性誤差理論對研究大氣和海洋科學中的非線性問題是非常必要的。Mu等[7-8]為了克服線性奇異向量方法的線性局限性，提出了條件非線性最優擾動（Conditional Nonlinear Optimal Perturbation, CNOP）方法。CNOP方法是線性奇異向量方法在非線性框架下自然的擴展。CNOP方法的建立對研究大氣和海洋科學的可預報性問題提供了有力的工具。

1 CNOP的定義

CNOP是在一定誤差范圍內，在預報時刻對預報結果不確定性產生最大影響的誤差。根據大氣和海洋科學可預報性的分類[9]，Mu等[7]建立了研究第一類可預報性問題的方法，與初始誤差有關的CNOP方法，記為CNOP-I方法。根據第二類可預報性問題，Mu等[8]建立了與初始誤差和模式誤差都有關的CNOP方法，其中與模式誤差有關的CNOP方法記為CNOP-P方法。具體地，假設一個非線性動力系統：

式中，X0表示狀態變量X的初始狀態。表示模式參數向量，m代表式（1）中參數的個數，F是一個非線性微分算子。在離散狀態下，式（1）在T時刻的解可以表示為：

這里MT(P)是在固定參數向量P時的非線性傳播算子，它將初始時刻的狀態X0“傳播”到T時刻的狀態X (T)。

用U0表示模式初值，u0表示模式的初始誤差，表示模式參數誤差，根據式（2），模式的解可寫為如下形式：

式中，度量了由聯合誤差導致狀態變量在T時刻偏離參考態U(T)的大小。

假設只考慮初始誤差u0，不考慮模式的參數誤差，即 =0，根據式（3），我們能獲得它所對應的T時刻的解U(T) +u(T) ，即有：

式中u(T)為該初始誤差的非線性發展。

假設只考慮模式參數誤差 ，不考慮模式的初始誤差，即u0=0，那么模式的解可寫為如下形式：

式中up(T)度量了由參數誤差導致狀態變量在T時刻偏離參考態U(T)的大小。

定義如下的非線性最優化問題：

式中，表示對初始誤差的約束，表示對參數誤差的約束，并且

顯然，最優化問題式（6）是一個有約束的最優化問題。對于給定的范數，式（6）中最優化問題的解（）為最優的初始誤差和參數誤差，表示在一定誤差約束條件下，在T時刻使得目標函數J偏離參考態程度最大。可以看到，當我們僅僅考慮初始誤差或者假設參數誤差 =0，最優化問題式（6）變為：

并且滿足式（8）的初始誤差u0δ就是CNOP-I，即為Mu等[7]定義的CNOP。如果在式（8）中用切線性模式的傳播算子代替MT，那么滿足式（8）式的初始誤差就是線性奇異向量（Linear Singular Vector, LSV），因此，CNOP方法是LSV方法在非線性框架下自然地推廣。

若在式（6）中我們不考慮初始誤差，即u0=0，可以得到：式中，是與參數有關的條件非線性最優擾動。類似于CNOP-I，稱這樣的最優參數誤差為CNOP-P。這樣可以明確地看到式（6）中的最優聯合模態是CNOP-I和CNOP-P的擴展（Mu等[8]）。

2 CNOP方法的應用

近些年，CNOP方法已經被應用于研究大氣和海洋科學中的一些物理問題[10]。例如，Mu等[11-12]探討了ENSO的可預報性和春季預報障礙問題（Spring Predictability Barrier, SPB）。研究結果表明CNOP型初始誤差有明顯的季節依賴性，在El Ni?o事件的不同發展階段，CNOP型初始誤差都會導致SPB。盡管LSV型初始誤差的發展也有季節依賴性，但在各個季節的增長率遠遠小于CNOP型初始誤差導致的預報誤差的增長率，這體現了CNOP和LSV間的差異[13-14]。Mu等[15]在熱鹽環流的非線性變化方面也進行了探索。研究結果指出通過對流過程，初始誤差的非線性作用導致了非對稱的響應。Jiang等[16]將CNOP方法也應用于集合預報的研究，結果表明由CNOP方法生成的集合樣本對中期集合預報的技巧高于由線性奇異向量方法。近些年，有一些學者將CNOP方法應用于臺風目標觀測的研究中。Mu等[17]研究結果表明用CNOP方法確定的敏感區進行目標觀測對臺風預報的改善程度大于用LSV確定的敏感區[18]。Terwiss-cha van Scheltinga等[19]用CNOP方法研究了正壓模式中雙環流的非線性變化，結果表明當雷諾數達到一定數值后，雙環流的非對稱且線性穩定的狀態將變為非線性不穩定性的。上述研究屬于第一類可預報性研究范疇。為了研究第二類可預報性問題，Duan等[20]用CNOP-P方法研究了ENSO事件的SPB現象。CNOP-P方法也在陸地生態系統變化研究中進行了應用。Sun等[21-23]將溫度和降水視為外強迫型參數，應用CNOP-P方法研究了溫度和降水變化情況下，中國區域陸地生態系統變化的對外強迫最大響應。結果表明在干旱和半干旱區域以及南方區域，陸地生態系統的響應較為明顯。

上述研究只是針對單獨考慮初始擾動或者模式參數擾動。近些年，一些研究同時考慮初始擾動和參數擾動。Wang等[24]用CNOP方法考察了模式參數誤差對黑潮路徑變異的影響，結果表明，盡管初始誤差導致的黑潮路徑變異的預報誤差大于模式誤差導致的預報誤差，但是模式誤差對黑潮路徑變異的影響不容忽視。Yu等[25]利用Zebiak-Cane模式分析了初始誤差和模式誤差在導致顯著春季預報障礙中所扮演的角色。數值結果表明在導致顯著春季預報障礙中，初始誤差比模式誤差更重要。Sun等[26]利用CNOP方法研究了草原生態系統的穩定性，結果表明在不同強度人類活動情況下，對草原生態系統模擬不確定性影響最大（即發生突變）的氣候變化類型是不相同的。這意味著氣候變化的模態對于草原生態系統的維持或退化是非常重要的，同時告訴我們應該合理地利用草原生態系統。上述這些研究說明了CNOP方法在大氣和海洋科學的非線性問題的研究中是一個有效的工具。

3 CNOP的數值求解

計算CNOP是一個求解有約束條件的非線性最優化問題的過程。如何計算這個非線性最優化問題是應用CNOP方法研究大氣和海洋科學的非線性物理問題的關鍵之處。我們以計算CNOP-I這個非線性最優化問題為例，計算CNOP-P與其是類似的。計算CNOP-I就是利用合適的最優化算法求解式（8）。事實上，已有很多較為成熟的最優化算法能夠求解公式（8）。序列二次規劃（Sequential Quadratic Programming, SQP，Barclay等[27]）算法﹑譜投影梯度（Spectral Projected Gradient, SPG，Birgin等[28]）算法和Limitedmemory Broyden-Fletcher-Goldfarb-Shanno（L-BFGS, Liu等[29]）方法等都是求解有約束的非線性最優化問題的最優化算法。如果能夠將目標函數﹑目標函數關于優化變量的梯度和約束條件準確地輸入到非線性最優化算法中，就可以通過計算得到CNOP-I。

目標函數﹑目標函數關于優化變量的梯度和約束條件需要使用者定義并且輸入到非線性最優化程序中（圖1）。首先，物理問題的約束條件可以較為容易地在約束優化算法中建立。有些約束優化算法也需要約束條件關于優化變量的梯度，例如SQP算法。這也可以較為容易地建立。其次，在優化算法中，需給出正確的目標函數值。式（7）中算子MT代表一個非線性模式從初始時刻到預報時刻T的傳播算子。為了得到目標函數，需要將非線性模式向前積分的信息輸入到優化算法中。即優化算法程序需要調用兩次模式的數值積分程序：一次是參考態（或基本態）MT(P) (U0)，此時并未考慮初始誤差；一次是考慮初始誤差時的非線性積分MT(P) (U0+u0)。將這兩次非線性積分的結果輸入到優化算法程序中，從而獲得目標函數。最后，目標函數關于優化變量的梯度是應用約束優化算法計算CNOP中非常重要的信息，它也是能否找到最優值的關鍵。定義法是求梯度最基本的方法。這種方法對理論模型是非常有效的。然而對于大氣和海洋科學中的海-氣耦合模式，由于該模式包含了較為復雜的動力框架和物理過程，用定義法求梯度的計算量非常大。對于大氣和海洋科學中的數值模式，獲得目標函數梯度的另一種技術是伴隨技術。在這里，我們不介紹伴隨技術。利用向前積分模式﹑其切線性模式和伴隨模式，可以獲得目標函數關于優化變量的梯度，將梯度信息輸入到約束優化算法中，即優化算法程序中需要調用向前積分模式﹑其切線性模式和伴隨模式三個程序。至此，上述三個信息輸入到約束優化算法中并且給出初始猜測值，算法就能正常運行并最終得到CNOP-I。

圖1 計算CNOP的流程圖Fig. 1 The flow chart to calculate CNOP

1）Burgers方程

Burgers方程如下所示：

它描述了狀態變量U隨時間的變化。在這里，γ=0.005m2/s，L=100m，T=30s。我們使用中心有限差分方法對式（10）進行離散化，得到了Burgers方程的數值積分模式，空間步長和時間步長分別是Δx=1.0m（空間格點總數目nx=101）、Δt=1.0s。然后在此基礎之上，寫出方程的切線性模式和伴隨模式（為了求目標函數關于優化變量的梯度）。使用Fortran進行編程，在附錄中我們給出了各模式的程序[30]。

2）切線性模式、伴隨模式以及梯度計算檢驗

在積分伴隨模式求目標函數關于優化變量的梯度之前，一般要對切線性模式、伴隨模式以及所得的梯度進行檢驗，從而保證實際應用中所得的結果與理論分析相符合。下面將分別給出這三部分的檢驗結果。需要說明的是：以下數值結果都是在變量和參數取雙精度的情況下得到的；各部分中擾動u0的取值都是u0(i)=Δ×i×0.01×(－1)i(i=1,…, nx－2)。Δ表示優化問題中的約束半徑。

① 切線性模式程序檢驗結果

根據切線性模式的定義，對所建立的切線性模式程序是否正確可以采用以下公式進行檢驗：

式中，M表示非線性模式，L表示M的切線性模式，表示取歐氏范數，u0是U0的擾動，0＜α＜1。某一基態下，切線性模式程序檢驗結果如表1所示。由此可知，當α趨于0時，R一致的逼近于1。這說明切線性模式的程序非常準確。

表1 切線性模式的檢驗結果Table1 The results of the tangent linear model

② 伴隨模式程序檢驗結果

根據伴隨算子的定義，伴隨模式程序的檢驗公式如下：

式中，L表示切線性模式，L*表示伴隨模式，u0代表擾動。

依照式（12），以u0為初值，積分切線性模式，將結果和其自身作內積，記作Valtlm。然后以Lu0為初值，積分伴隨模式，將結果和u0作內積，記作Valadj。最后在機器精度的標準下，檢驗Valtlm=Valadj是否成立。某一基態下，伴隨模式檢驗結果如表2所示。

表2 伴隨模式的檢驗結果Table2 The results of the adjoint model

由表2可知，Valtlm和Valadj相等的有效數字達到15位，非常接近機器精度（16位有效數字），這說明從切線性模式程序出發建立的伴隨模式程序是準確的。

③ 梯度計算檢驗結果

根據定義，梯度計算的檢驗公式如下：

式中，J(u0)是目標函數，是由伴隨程序計算的J在擾動u0處的梯度表示取歐氏范數，0＜α＜1。梯度的檢驗結果如表3所示。由表可知，當α趨于0時，R在1附近呈現拋物型結構，這說明伴隨模式可用于求解優化問題中代價函數的梯度。

表3 梯度計算的檢驗結果Table3 The results of gradients

3）計算結果

圖2是當Burgers方程中初始條件為U0時，狀態變量U隨時間的發展，也是我們求解CNOP時的基態，也可以稱其為參考態，即式（7）中MT(P) (U0)。

圖2 初始條件為U0時狀態變量U的非線性發展Fig. 2 The nonlinear evolution of the state variableUwhen the initial condition isU0

在求解目標函數J(u0)最小值和最小值點（CNOP）時，我們使用的優化算法是SPG2，該算法需要使用目標函數的梯度信息。下面分別使用定義法和伴隨方法獲取目標函數的梯度，進而求得CNOP，并且比較兩種方法得到的CNOP有何差異。

如圖3所示，兩種方法所得到的CNOP空間分布極其相似。圖4是兩個CNOP差異的空間分布。由上可知，基于伴隨方法得到的CNOP和基于定義方法得到的CNOP幾乎是一樣的，而且兩者隨時間的發展也幾乎一樣，如圖5所示。綜上可知，兩種方法得到的CNOP幾乎一樣。然而基于伴隨求CNOP的方法所使用的CPU時間遠遠小于基于定義求CNOP的方法所使用的CPU時間（圖5）。因而，基于伴隨求CNOP的方法準確而且效率高。圖6給出了CNOP非線性發展的時空分布。由此看到在時刻T，CNOP類型初始擾動有最大非線性發展。

圖3 CNOP的空間分布（a）基于伴隨方法；（b）基于定義方法Fig. 3 The spatial distribution of CNOP（a）The adjoint method；（b）The definition method

圖4 基于伴隨方法和定義方法計算得到的CNOP的差異Fig. 4 The difference of CNOP obtained by the adjoint method and the definition method

圖5 基于伴隨方法和定義方法的CNOP的非線性發展Fig. 5 The nonlinear evolution of CNOP obtained by the adjoint method and the definition method

圖6 CNOP非線性演變的時空分布特征Fig. 6 The spatial and temporal difference of the nonlinear evolution of CNOP

4 結論和討論

天氣和氣候的可預報性問題是目前科學研究的熱點問題之一。傳統的LSV方法在研究天氣和氣候的可預報性問題上具有局限性，因為此方法是假設天氣和氣候事件的演變過程是線性的。為了擴展LSV方法的局限性，Mu等[7-8]建立了CNOP方法。本文不僅介紹了CNOP的定義和拓展CNOP等定義，并且詳細地介紹了CNOP方法在ENSO、黑潮和阻塞可預報性，以及熱鹽環流和草原生態系統穩定性研究中的應用。

特別地，為了讓更多的學者能夠應用CNOP方法并將其應用于更多的研究領域，本文利用一個理論的Burgers方程，詳細地介紹了如何建立此方法的積分模式、切線性模式和伴隨模式[30]，通過介紹計算CNOP的流程圖，嘗試讓學者了解如何利用非線性最優化算法計算CNOP。本文只是介紹了與初始擾動有關的CNOP的計算，如果計算CNOP-P，可以參考Mu等[8]的研究。從本文的介紹可以看出，計算CNOP需要積分模式對應的切線性模式和伴隨模式。然而，目前有些海-氣耦合模式缺乏相應的切線性模式和伴隨模式。為了能使CNOP方法得到更廣泛的應用，一些學者[31-32]利用集合投影方法近似計算梯度信息，從而避免了需要積分模式對應的切線性模式和伴隨模式而獲得梯度信息。這為CNOP方法的應用提供了一種可行的途徑。

[1]Leith C E. Theoretical skill of Monte Carlo forecast. Mon Wea Rev, 1974, 102: 409-418.

[2]Houtekamer P L, Lefaivre L, Derome J, et al. A system simulation approach to ensemble prediction. Mon Wea Rev, 1996, 124: 1225-1242.

[3]Toth Z, Kalnay E. Ensemble forecasting at NCEP: the breeding method. Mon Wea Rev, 1997, 125: 3297-3318.

[4]Buizza R, Palmer T. Potential forecast skill of ensemble prediction, and spread and skill distributions of the ECMWFensemble prediction system. Mon Wea Rev, 1997, 125: 99-119.

[5]Hamill T M, Snyder C, Morss R E. A comparison of probabilistic forecasts from bred, singular-vector, and perturbed observation ensembles. Mon Wea Rev, 2000, 128: 1835-1851.

[6]Garcia-Moya J A, Callado A, Escriba P, et al. Predictability of short-range forecasting: a multimodel approach. Tellus A, 2011, 63: 550-563.

[7]Mu M, Duan W S, Wang B. Conditional nonlinear optimal perturbation and its applications. Nonlin Processes Geophys , 2003, 10: 493-501.

[8]Mu M, Duan W S, Wang Q, et al. An extension of conditional nonlinear optimal perturbation approach and its applications. Nonlin Processes Geophys, 2010, 17: 21-220. doi:10.5194/ npg-17-211-2010.

[9]Mu M, Duan W S, Wang J C. The predictability problems in numerical weather and climate prediction. Adv Atmos Sci, 2002, 19(2): 191-204.

[10]穆穆,段晚鎖.2013.條件非線性最優擾動在可預報性問題研究中的應用.大氣科學,37(2):281-296. doi:10.3878/ j.issn.1006-9895.2012.12319.

[11]Mu M, Duan W S, Wang B. Season-dependent dynamics of nonlinear optimal error growth and El Ni?o-Southern Oscillation predictability in a theoretical model. J Geophys Res, 2007, 112, D10113. doi:10.1029/2005JD006981.

[12]Mu M, Xu H, Duan W S. A kind of initial errors related to “spring predictability barrier” for El Ni?o events in Zebiak-Cane model. Geophys Res Lett, 2007, 34, L03709. doi:10.1029/2006GL027412.

[13]Duan W S, Xu H, Mu M. Decisive role of nonlinear temperature advection in El Ni?o- and La Ni?a- amplitude asymmetry. J Geophys Res, 2008, 113, C01014. doi:10.1029/2006JC003974.

[14]Mu M, Duan W S, Wang J F. Nonlinear optimization problems in atmospheric and oceanic sciences. East-west J Math, Thailand, 2002, (S): 155-164.

[15]Mu M, Sun L, Dijikstra H A. The sensitivity and stability of the ocean’s thermohaline circulation to fi nite amplitude perturbations. J Phys Oceanogr, 2004, 34: 2305-2315.

[16]Jiang Z N, Mu M. A comparisons study of the methods of conditional nonlinear optimal perturbations and singular vectors in ensemble prediction. Adv Atmos Sci, 2009, 26(3): 465-470. doi: 10.1007/s00376-009-0465-6.

[17]Mu M, Zhou F F, Wang H L. A method for identifying the sensitive areas in targeted observations for tropical cyclone prediction: conditional nonlinear optimal perturbation. Mon Wea Rev, 2009, 137: 1623-1639.

[18]Qin X, Mu M. Influence of conditional nonlinear optimal perturbations sensitivity on typhoon track forecasts. Quart J Roy Meteor Soc, 2011, 138: 185-197.

[19]Terwisscha van S A D, Dijkstra H A. Conditional nonlinear optimal perturbations of the double-gyre ocean circulation. Nonlin Processes Geophys, 2008, 15: 727-734. doi:10.5194/npg-15-727-2008.

[20]Duan W S, Zhang R. Is model parameter error related to a signi fi cant spring predictability barrier for El Nino events? Results from a theoretical model. Adv Atmos Sci, 2010, 27(5): 1003-1013.

[21]Sun G D, Mu M. Responses of soil carbon variation to climate variability in China using the LPJ model. Theoretical and Applied Climatology, 2012, 110:143-153. doi:10.1007/s00704-012-0619-9. [22]Sun G D, Mu M. Understanding variations and seasonal characteristics of net primary production under two types of climate change scenarios in China using the LPJ model. Climatic Change, 2013, 120:755-769.doi: 10.1007/s10584-013-0833-1.

[23]Sun G D, Mu M. The analyses of the net primary production due to regional and seasonal temperature di ff erences in eastern China using the LPJ model. Ecological Modelling, 2014, 289: 66-76. doi: 10.1016/j.ecolmodel.2014.06.021.

[24]Wang Q, Mu M, Dijkstra H A. Application of the conditional nonlinear optimal perturbation method to the predictability study of the Kuroshio large meander. Adv Atmos Sci, 2012, 29(1): 118-134. doi: 10.1007/s00376-011-0199-0.

[25]Yu Y S, Mu M, Duan W S. Does Model Parameter Error Cause a Signi fi cant “Spring Predictability Barrier” for El Ni?o Events in the Zebiak–Cane Model? J Climate, 2012, 25: 1263-1277. doi: http://dx.doi.org/10.1175/2011JCLI4022.1

[26]Sun G D, Mu M. Inducing unstable grassland equilibrium states due to nonlinear optimal patterns of initial and parameter perturbations: theoretical models. Adv Atmos Sci, 2012, 29(1): 79-90. doi: 10.1007/s00376-011-0226-1.

[27]Barclay A, Gill P E, Rosen J B. SQP methods and their application to numerical optimal control. Numerical Analysis Report 97-3, Department of Mathematics, University of California, San Diego, La Jolla, CA, 1997.

[28]Birgin E G, Martinez J M, Raydan M. Nonmonotone spectral projected gradient methods on convex sets. SIAM J Optimization, 2000, 10(4): 1196-1211.

[29]Liu D C, Nocedal J. On the limited memory BFGS method for large scale minimization. Math Prog, 1989, 45: 503-528.

[30]Kalnay E. Atmospheric Modeling, Data Assimilation, and Predictability. Cambridge: Cambridge University Press, 2003.

[31]Wang B, Tan X W. Conditional nonlinear optimal perturbations: adjoint-free calculation method and preliminary test. Mon Wea Rev, 2010, 138: 1043-1049

[32]Chen L, Duan W S, Xu H. A SVD-based ensemble projection algorithm for calculating conditional nonlinear optimal perturbation. Sci China: Earth Sci,2014,58 :385-394. doi: 10.1007/ s11430-014-4991-4.

附錄1 Burgers方程的數值積分模式、切線性模式和伴隨模式

1）數值積分模式

SUBROUTINE BURGER(nx,n,ui,ub,u)

implicit none

integer nx !number of grid points

integer n !number of time steps

double precision ub(nx,n) !basic states(u)

double precision u(nx,n) !model solutions

double precision ui(nx) !initial conditions

double precision dx !space increment

double precision dt !time increment

double precision gamma ! diffusion coeff i cient

double precision c0,c1

integer i,j

common /com_param/ c0,c1

gamma=0.005d0

dx=0.01d0

dt=0.001d0

c0=dt/dx

c1=gamma*dt/(dx**2)

!set the initial conditions:

do i=1,nx

u(i,1)=ui(i)

end do

!set the boundary conditions:

do j=1,n

u(1,j)=0.

u(nx,j)=0.

end do

!integerate the model numerically:-----use central fi nite difference method on spatial and temporal directions

do i=2,nx-1

u(i,2)=u(i,1)-0.25*c0*(u(i+1,1)*u(i+1,1)-u(i-1,1)*u(i-1,1))+c1*(u(i+1,1)+u(i-1,1)-2*u(i,1))

end do

do j=3,n

do i=2,nx-1

u(i,j)=u(i,j-2)-0.5*c0*(u(i+1,j-1)*u(i+1,j-1)-u(i-1,j-1)*u(i-1,j-1))+2*c1*(u(i+1,j-1)+u(i-1,j-1)-2*u(i,j-1))

end do

end do

!save the nonlinear solutions to the basic fi elds:

do j=1,n

do i=1,nx

ub(i,j)=u(i,j)

end do

end do return

END SUBROUTINE

2）切線性模式

SUBROUTINE BURGER_TLM(nx,n,ui,ub,u)

implicit none

integer nx,n

double precision ub(nx,n) !basic state

double precision u(nx,n) !TLM solutions

double precision ui(nx) !initial conditions

common /com_param/ c0,c1

double precision c0,c1

integer i,j

! set the initial conditions

do i=1,nx

u(i,1)=ui(i)

end do

! set the bundary conditions:

do j=1,n

u(1,j)=0.

u(nx,j)=0.

end do

! integer the model numerically:

do i=2,nx-1

u(i,2)=u(i,1)-0.5*c0*(ub(i+1,1)*u(i+1,1)-ub(i-1,1)*u(i-1,1))+c1*(u(i+1,1)+u(i-1,1)-2*u(i,1))

end do

do j=3,n

do i=2,nx-1

u(i,j)=u(i,j-2)-c0*(ub(i+1,j-1)*u(i+1,j-1)-ub(i-1,j-1)*u(i-1,j-1))+2*c1*(u(i+1,j-1)+u(i-1,j-1)-2*u(i,j-1))

end do

enddo

! set the fi nal value of u in ui

do i=1,nx

ui(i)=u(i,n)

end do

return

END SUBROUTINE

3) 伴隨模式

SUBROUTINE BURGER_ADJ(nx,n,ui,ub,u)

implicit none

integer nx,n

double precision ub(nx,n) !basic states

double precision u(nx,n) !modle solutions

double precision ui(nx) !initial conditions

double precision dt !time increment

double precision dx !space increment

double precision c0,c1

integer i,j

common /com_param/ c0,c1

!initialize adjoint variables:

do j=1,n

do i=1,nx

u(i,j)=0.

end do

end do

!set the fi nal conditions:

do i=1,nx

u(i,n)=ui(i)

ui(i)=0.

end do

do j=n,3,-1

do i=nx-1,2,-1

u(i,j-2)=u(i,j-2)+u(i,j)

u(i-1,j-1)=u(i-1,j-1)+(2*c1+c0*ub(i-1,j-1))*u(i,j)

u(i,j-1)=u(i,j-1)-4*c1*u(i,j)

u(i+1,j-1)=u(i+1,j-1)+(2*c1-c0*ub(i+1,j-1))*u(i,j)

u(i,j)=0

end do

end do

do i=nx-1,2,-1

u(i-1,1)=u(i-1,1)+(c1+0.5*c0*ub(i-1,1))*u(i,2)

u(i,1)=u(i,1)+(1-2*c1)*u(i,2)

u(i+1,1)=u(i+1,1)+(c1-0.5*c0*ub(i+1,1))*u(i,2) u(i,2)=0

end do

!set the boundary conditions:

do i=1,n

u(1,i)=0.0

u(nx,i)=0.

end do

!set the fi nal value of u in ui:----Gradent

do i=1,nx

ui(i)=ui(i)+u(i,1)

end do

return

END SUBROUTINE

附錄2 譜投影梯度（SPG2）算法簡介

譜投影梯度（Spectral Projected Gradient Method，SPG2）是一種非單調梯度投影算法，用于解決大規模的、可行域是凸集的優化問題（Birgin 等）[A1]。這種優化問題具有如下形式：

式中，Ω是Rn中的閉凸集，而且目標函數f在包含Ω的開集上具有連續偏導數。

需要說明的是，對于任意的，定義是與給定范數有關的在閉凸集Ω上的投影。g(x)表示目標函數的梯度，即。該算法以任意的開始。整數m表示算法中存儲的目標函數值的最大個數；αmin是一個小參數，并且大于0；αmax是一個大參數，αmax＞αmin；，代表一個足夠小的參數；參數σ1，σ2滿足：0＜σ1＜σ2＜1；表示內積。SPG2算法具體如下：

Step1 計算譜投影梯度。檢驗算法終止準則是否滿足。如果滿足，算法結束，xk即為最小值點；否則執行Step2。

Step4 迭代次數k=k+1。進行Step1。

在SPG2優化算法中，收斂準則是：或者，其中ε1，ε2都是很小的數。如果收斂準則滿足，則算法終止。此外，如果迭代次數k超過了給定的最大迭代次數maxit，或者目標函數值的計算次數超過了給定的最大計算次數maxfun，則算法同樣會終止。此外，需要注意的是這個算法中目標函數值和目標函數的梯度需要由使用者給出。

參考文獻

[A1] Birgin E G, Martinez J M, Raydan M. SPG: software for convexconstrained optimization. ACM Trans Math Softw, 2001, 27(3)：340-349.

Conditional Nonlinear Optimal Perturbation: Introduction and Numerical Computation

Sun Guodong1,4, Mu Mu2, Duan Wansuo1, Wang Qiang3, Peng Fei1,4
(1 State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics (LASG), Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029 2 Institute of Atmospheric Sciences, Fudan University, Shanghai 200433 3 Key Laboratory of Ocean Circulation and Wave, Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071 4 University of Chinese Academy of Sciences, Beijing 100049)

This paper introduces the def i nition of conditional nonlinear optimal perturbation (CNOP), and the applications of the CNOP in atmosphere and ocean studies. The CNOP approach is expanded as that related to initial perturbation (CNOP-I), related to parameter perturbation (CNOP-P), and the combined both of CNOP-I and CNOP-P, according to the different perturbation types. The CNOP-I approach has been applied to the predictability studies of ENSO events, Kuroshio path anomalies, blocking, nonlinear stabilities of thermohaline circulation and grassland ecosystem. The CNOP-I has been further employed to explore the target observation of typhoon. The sensitive region could be identif i ed by using the CNOP-I approach. The forecast skill may be improved by adding more adaptive observations in the sensitive region. The CNOP-P approach has been applied also to Kuroshio path anomalies, nonlinear stabilities of thermohaline circulation and grassland ecosystem. Here, we carried out a numerical simulation how to obtain the CNOP with the Burgers equation through building the tangent linear model and adjoint model. The result shows that the CNOP can be calculated by using the Burgers equation, the tangent linear model and the adjoint model with nonlinear optimization algorithm; It supplies a guide to a beginner to learn the CNOP and a reference for employing the CNOP to other applicable subjects.

conditional nonlinear optimal perturbation (CNOP), predictability, adaptive observations

10.3969/j.issn.2095-1973.2016.06.001

2014年10月23日；

2015年3月31日

孫國棟（1980 —），Email: sungd@mail.iap.ac.cn

資助信息：國家自然科學基金（41375111，91437111）