孫國(guó)棟穆穆段晚鎖王強(qiáng)彭飛

（1 中國(guó)科學(xué)院大氣物理研究所，大氣科學(xué)和地球流體力學(xué)數(shù)值模擬國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029；2 復(fù)旦大學(xué)大氣科學(xué)研究院，上海 200433；3 中國(guó)科學(xué)院海洋研究所，海洋環(huán)流與波動(dòng)重點(diǎn)實(shí)驗(yàn)室，青島 266071；4 中國(guó)科學(xué)院大學(xué)，北京 100049）

條件非線性最優(yōu)擾動(dòng)（CNOP）：簡(jiǎn)介與數(shù)值求解

孫國(guó)棟1,4穆穆2段晚鎖1王強(qiáng)3彭飛1,4

（1 中國(guó)科學(xué)院大氣物理研究所，大氣科學(xué)和地球流體力學(xué)數(shù)值模擬國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029；2 復(fù)旦大學(xué)大氣科學(xué)研究院，上海 200433；3 中國(guó)科學(xué)院海洋研究所，海洋環(huán)流與波動(dòng)重點(diǎn)實(shí)驗(yàn)室，青島 266071；4 中國(guó)科學(xué)院大學(xué)，北京 100049）

介紹了條件非線性最優(yōu)擾動(dòng)（Conditional Nonlinear Optimal Perturbation，CNOP）的定義及其在大氣和海洋等可預(yù)報(bào)性研究中的應(yīng)用。根據(jù)研究對(duì)象不同，CNOP分為與初始擾動(dòng)有關(guān)的CNOP（CNOP-I）方法、與模式參數(shù)擾動(dòng)有關(guān)的CNOP（CNOP-P）方法和同時(shí)考慮初始擾動(dòng)和模式參數(shù)擾動(dòng)的CNOP方法。目前，CNOP-I方法已經(jīng)應(yīng)用于ENSO、黑潮和阻塞可預(yù)報(bào)性以及熱鹽環(huán)流和草原生態(tài)系統(tǒng)穩(wěn)定性的研究。此外，CNOP-I方法也被應(yīng)用于探討臺(tái)風(fēng)目標(biāo)觀測(cè)的研究，利用CNOP-I方法能夠識(shí)別出臺(tái)風(fēng)預(yù)報(bào)的初值敏感區(qū)，通過(guò)觀測(cè)系統(tǒng)模擬試驗(yàn)表明在初值敏感區(qū)增加觀測(cè)能夠有效改進(jìn)臺(tái)風(fēng)的預(yù)報(bào)技巧。CNOP-P方法也在ENSO和黑潮可預(yù)報(bào)性以及熱鹽環(huán)流和草原生態(tài)系統(tǒng)穩(wěn)定性研究中得到了應(yīng)用。為了將CNOP方法應(yīng)用于更多的領(lǐng)域，本文利用一個(gè)簡(jiǎn)單的Burgers方程，介紹了如何通過(guò)建立Burgers方程的切線性模式和伴隨模式，從而利用非線性最優(yōu)化算法計(jì)算獲得CNOP。這一數(shù)值試驗(yàn)為將CNOP方法應(yīng)用于更多的領(lǐng)域提供了借鑒。

條件非線性最優(yōu)擾動(dòng)方法（CNOP），可預(yù)報(bào)性，目標(biāo)觀測(cè)

0 引言

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

1 CNOP的定義

CNOP是在一定誤差范圍內(nèi)，在預(yù)報(bào)時(shí)刻對(duì)預(yù)報(bào)結(jié)果不確定性產(chǎn)生最大影響的誤差。根據(jù)大氣和海洋科學(xué)可預(yù)報(bào)性的分類[9]，Mu等[7]建立了研究第一類可預(yù)報(bào)性問(wèn)題的方法，與初始誤差有關(guān)的CNOP方法，記為CNOP-I方法。根據(jù)第二類可預(yù)報(bào)性問(wèn)題，Mu等[8]建立了與初始誤差和模式誤差都有關(guān)的CNOP方法，其中與模式誤差有關(guān)的CNOP方法記為CNOP-P方法。具體地，假設(shè)一個(gè)非線性動(dòng)力系統(tǒng)：

式中，X0表示狀態(tài)變量X的初始狀態(tài)。表示模式參數(shù)向量，m代表式（1）中參數(shù)的個(gè)數(shù)，F(xiàn)是一個(gè)非線性微分算子。在離散狀態(tài)下，式（1）在T時(shí)刻的解可以表示為：

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

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

式中，度量了由聯(lián)合誤差導(dǎo)致?tīng)顟B(tài)變量在T時(shí)刻偏離參考態(tài)U(T)的大小。

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

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

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

式中up(T)度量了由參數(shù)誤差導(dǎo)致?tīng)顟B(tài)變量在T時(shí)刻偏離參考態(tài)U(T)的大小。

定義如下的非線性最優(yōu)化問(wèn)題：

式中，表示對(duì)初始誤差的約束，表示對(duì)參數(shù)誤差的約束，并且

顯然，最優(yōu)化問(wèn)題式（6）是一個(gè)有約束的最優(yōu)化問(wèn)題。對(duì)于給定的范數(shù)，式（6）中最優(yōu)化問(wèn)題的解（）為最優(yōu)的初始誤差和參數(shù)誤差，表示在一定誤差約束條件下，在T時(shí)刻使得目標(biāo)函數(shù)J偏離參考態(tài)程度最大?？梢钥吹?，當(dāng)我們僅僅考慮初始誤差或者假設(shè)參數(shù)誤差 =0，最優(yōu)化問(wèn)題式（6）變?yōu)椋?/p>

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

若在式（6）中我們不考慮初始誤差，即u0=0，可以得到：式中，是與參數(shù)有關(guān)的條件非線性最優(yōu)擾動(dòng)。類似于CNOP-I，稱這樣的最優(yōu)參數(shù)誤差為CNOP-P。這樣可以明確地看到式（6）中的最優(yōu)聯(lián)合模態(tài)是CNOP-I和CNOP-P的擴(kuò)展（Mu等[8]）。

2 CNOP方法的應(yīng)用

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

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

3 CNOP的數(shù)值求解

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

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

圖1 計(jì)算CNOP的流程圖Fig. 1 The flow chart to calculate CNOP

1）Burgers方程

Burgers方程如下所示：

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

2）切線性模式、伴隨模式以及梯度計(jì)算檢驗(yàn)

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

① 切線性模式程序檢驗(yàn)結(jié)果

根據(jù)切線性模式的定義，對(duì)所建立的切線性模式程序是否正確可以采用以下公式進(jìn)行檢驗(yàn)：

式中，M表示非線性模式，L表示M的切線性模式，表示取歐氏范數(shù)，u0是U0的擾動(dòng)，0＜α＜1。某一基態(tài)下，切線性模式程序檢驗(yàn)結(jié)果如表1所示。由此可知，當(dāng)α趨于0時(shí)，R一致的逼近于1。這說(shuō)明切線性模式的程序非常準(zhǔn)確。

表1 切線性模式的檢驗(yàn)結(jié)果Table1 The results of the tangent linear model

② 伴隨模式程序檢驗(yàn)結(jié)果

根據(jù)伴隨算子的定義，伴隨模式程序的檢驗(yàn)公式如下：

式中，L表示切線性模式，L*表示伴隨模式，u0代表擾動(dòng)。

依照式（12），以u(píng)0為初值，積分切線性模式，將結(jié)果和其自身作內(nèi)積，記作Valtlm。然后以Lu0為初值，積分伴隨模式，將結(jié)果和u0作內(nèi)積，記作Valadj。最后在機(jī)器精度的標(biāo)準(zhǔn)下，檢驗(yàn)Valtlm=Valadj是否成立。某一基態(tài)下，伴隨模式檢驗(yàn)結(jié)果如表2所示。

表2 伴隨模式的檢驗(yàn)結(jié)果Table2 The results of the adjoint model

由表2可知，Valtlm和Valadj相等的有效數(shù)字達(dá)到15位，非常接近機(jī)器精度（16位有效數(shù)字），這說(shuō)明從切線性模式程序出發(fā)建立的伴隨模式程序是準(zhǔn)確的。

③ 梯度計(jì)算檢驗(yàn)結(jié)果

根據(jù)定義，梯度計(jì)算的檢驗(yàn)公式如下：

式中，J(u0)是目標(biāo)函數(shù)，是由伴隨程序計(jì)算的J在擾動(dòng)u0處的梯度表示取歐氏范數(shù)，0＜α＜1。梯度的檢驗(yàn)結(jié)果如表3所示。由表可知，當(dāng)α趨于0時(shí)，R在1附近呈現(xiàn)拋物型結(jié)構(gòu)，這說(shuō)明伴隨模式可用于求解優(yōu)化問(wèn)題中代價(jià)函數(shù)的梯度。

表3 梯度計(jì)算的檢驗(yàn)結(jié)果Table3 The results of gradients

3）計(jì)算結(jié)果

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

圖2 初始條件為U0時(shí)狀態(tài)變量U的非線性發(fā)展Fig. 2 The nonlinear evolution of the state variableUwhen the initial condition isU0

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

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

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

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

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

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

4 結(jié)論和討論

天氣和氣候的可預(yù)報(bào)性問(wèn)題是目前科學(xué)研究的熱點(diǎn)問(wèn)題之一。傳統(tǒng)的LSV方法在研究天氣和氣候的可預(yù)報(bào)性問(wèn)題上具有局限性，因?yàn)榇朔椒ㄊ羌僭O(shè)天氣和氣候事件的演變過(guò)程是線性的。為了擴(kuò)展LSV方法的局限性，Mu等[7-8]建立了CNOP方法。本文不僅介紹了CNOP的定義和拓展CNOP等定義，并且詳細(xì)地介紹了CNOP方法在ENSO、黑潮和阻塞可預(yù)報(bào)性，以及熱鹽環(huán)流和草原生態(tài)系統(tǒng)穩(wěn)定性研究中的應(yīng)用。

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

[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.條件非線性最優(yōu)擾動(dòng)在可預(yù)報(bào)性問(wèn)題研究中的應(yīng)用.大氣科學(xué),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方程的數(shù)值積分模式、切線性模式和伴隨模式

1）數(shù)值積分模式

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）算法簡(jiǎn)介

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

式中，Ω是Rn中的閉凸集，而且目標(biāo)函數(shù)f在包含Ω的開(kāi)集上具有連續(xù)偏導(dǎo)數(shù)。

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

Step1 計(jì)算譜投影梯度。檢驗(yàn)算法終止準(zhǔn)則是否滿足。如果滿足，算法結(jié)束，xk即為最小值點(diǎn)；否則執(zhí)行Step2。

Step4 迭代次數(shù)k=k+1。進(jìn)行Step1。

在SPG2優(yōu)化算法中，收斂準(zhǔn)則是：或者，其中ε1，ε2都是很小的數(shù)。如果收斂準(zhǔn)則滿足，則算法終止。此外，如果迭代次數(shù)k超過(guò)了給定的最大迭代次數(shù)maxit，或者目標(biāo)函數(shù)值的計(jì)算次數(shù)超過(guò)了給定的最大計(jì)算次數(shù)maxfun，則算法同樣會(huì)終止。此外，需要注意的是這個(gè)算法中目標(biāo)函數(shù)值和目標(biāo)函數(shù)的梯度需要由使用者給出。

參考文獻(xiàn)

[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日

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

資助信息：國(guó)家自然科學(xué)基金（41375111，91437111）