耿煥同,孫義杰,張建,閔錦忠

(南京信息工程大學1.氣象災害省部共建教育部重點實驗室;2.計算機與軟件學院,江蘇南京 210044)

0 引言

氣象學中的天氣預報是在一定的模式大氣條件下,根據(jù)初值條件,用數(shù)值方法求解模式方程組的過程。造成預報誤差的主要原因是初始條件的不準確和預報模式的不準確,因此模式的不斷完善和高分辨觀測資料的不斷利用是推動天氣預報發(fā)展的主要途徑[1-2]。在實際預報業(yè)務中,往往先利用已有物理規(guī)律和觀測資料,對模式參數(shù)調整和初始場修正屬微分方程組求解中反問題研究范疇;真正作數(shù)值天氣預報時,實為正問題的求解[3]。

針對數(shù)值天氣預報中模式參數(shù)反演問題,傳統(tǒng)方法是通過經(jīng)驗或多次實驗去調整參數(shù),因此參數(shù)選擇必存在一定的隨意性和片面性,難以保證所作的調整是有效的;一旦某些參數(shù)選擇不當,往往會使預報準確率明顯降低[4]。其他傳統(tǒng)的一些求解反問題的方法存在限制,如增強型拉格朗日法對問題本身有比較苛刻的要求(如連續(xù)、可微)[5]。鑒此尋找一種新的適合于氣象上的微分方程參數(shù)識別反問題求解方法顯得十分必要。

隨著進化算法理論研究和工程應用的不斷深入和發(fā)展,進化算法(EA)已成為主要的優(yōu)化求解方法之一[6-7]。進化算法是一種模擬生物進化規(guī)律的搜索和優(yōu)化計算方法,如遺傳算法、進化策略等;此類算法具有很強的通用性、魯棒性,并且便于并行處理等特點,已被成功應用于機器學習、模式識別、優(yōu)化控制等領域中[8-9]。同時,在其他工程領域,應用進化算法進行反問題求解已取得許多研究成果[10-11]。

本文利用進化策略在微分方程組反問題求解中的優(yōu)勢,結合氣象問題,提出了一種基于進化策略的氣象學反問題求解算法,用于預報模式中參數(shù)識別及初始條件優(yōu)化兩類反問題的求解中。當利用微分方程作為預報模式進行實際預報時,利用已有觀測數(shù)據(jù)對原模式的參數(shù)(或初始場)進行反演以便改善預報水平;最后運用該方法對一維擴散方程和Lorenz-96簡單預報模式中反問題進行模擬實驗。

1 問題描述

1.1 微分方程反問題

在自然科學和工程應用中,許多模型都可以用微分方程來描述。微分方程的形成和發(fā)展是和力學、天文學、物理學,以及其他科學技術的發(fā)展密切相關的。通常微分方程的求解是尋找滿足定解問題初、邊值條件的微分方程的解,又稱為正問題。其特點是由“原因”推出“結果”,是自然探索的一種傳統(tǒng)途徑。微分方程的反問題是指由“結果”推出“原因”。即已知微分方程的解,去求方程中的未知部分。反問題是從各個領域、學科的實際需求中提出的,反問題研究是一門交叉性學科。例如,工程應用中的識別、反演問題就屬于反問題的范疇。

1.2 微分方程反問題的分類

設研究的微分方程一般形式如下:

微分方程Ku(x,t)=f(x,t),x∈Ω,t∈(0,∞)。

初始條件Iu(x,0)=φ(x),x∈Ω。

邊界條件Bu(x,t)=Φ(x,t),x∈aΩ。

附加條件Au(x,t)=Ψ(x,t);x∈aΩ。

其中:u(x,t)為偏微分方程的解;f(x,t)為右端項;φ(x)、Φ(x,t)、Ψ(x,t)分別為初始條件、邊界條件和附加條件。K、I、B、A分別為微分算子、初始算子、邊界算子和附加算子。當上述這些已知量中有一個變?yōu)槲粗獣r,就構成了偏微分方程的反問題。

1)當微分算子K為未知時,稱為算子識別反問題。通常情況下算子結構已知,未知的是算子中的參數(shù)。這類問題稱為參數(shù)反演、參數(shù)識別問題。在自然科學和實際應用中,這類問題最為常見。

2)當右端項f(x,t)為未知時,稱為源的反問題。

3)當初始條件φ(x)未知時,稱為初始條件反問題。在附加條件中往往給出系統(tǒng)在某一時刻的狀態(tài),即在后某一時刻去確定初始的狀態(tài),所以也稱為逆時間過程的反問題。

4)當邊界條件Φ(x,t)未知時,工程上稱為邊界控制反問題。

5)當區(qū)域的邊界aΩ為未知時,稱之類問題為幾何反問題。工程中的某些定向設計問題常歸為這類問題。一維幾何反問題是確定端點的位置,二維、三維的反問題則是確定曲線的形狀。

2 應用進化策略求解微分方程反問題

2.1 微分方程中的反問題求解到函數(shù)優(yōu)化的轉化

以基于以上定義的偏微分方程(除若干參數(shù)待調整外)預報系統(tǒng)為例,只要輸入初始場及模式參數(shù)值,模式即可輸出未來某一時刻的預報,以向量Bm表示t=tm時刻的初始場,實向量K表示模式參數(shù),向量Om表示t=tm+τ時刻模式的輸出,則預報過程可表示為

若M模式已知,但模式中的參數(shù)K未知,即為參數(shù)反演問題。利用觀測場O(或其中部分元素)的觀測值Oob,在參數(shù)K的容許區(qū)間內找到最佳近似Kbest,使模式輸出與相應的觀測值之差最小。

其中:‖‖表示兩個向量間的距離函數(shù),將‖M(K,B)-Oob‖作為目標函數(shù),即將參數(shù)反演問題轉化為函數(shù)優(yōu)化問題,根據(jù)提供的觀測值對模式中的參數(shù)進行調整,力使‖M(K,B)-Oob‖取得最小值,從而得到參數(shù)的最優(yōu)估計值。

若初始場B未知,即為初始場反問題。利用觀測場O(或其中部分元素)的觀測值Oob在觀測變量參數(shù)B的容許區(qū)間內找到最佳近似初始場(或分析場)Ba,使模式輸出與相應的觀測誤差最小。

將‖M(K,B)-Oob‖作為目標函數(shù),即將初始場反問題轉化為函數(shù)優(yōu)化問題,根據(jù)觀測值調整參數(shù)B,使‖M(K,B)-Oob‖取得最小值時,就得到初始場的最優(yōu)估計值。若觀測變量間需滿足協(xié)調條件,則可轉化為帶約束的優(yōu)化問題,本文暫不作討論。

綜上所述,令目標函數(shù)f=inf‖M(K,B)-Oob‖,將反問題的求解等價轉化為函數(shù)優(yōu)化問題。即根據(jù)所提供的觀測資料采用進化策略優(yōu)化方法,調整模式參數(shù)(或初始場Ba)使‖M(K,B)-Oob‖取得最小值,以得到模式參數(shù)(或初始場Bm)的最優(yōu)估計值。

2.2 進化策略

經(jīng)過從反問題求解到函數(shù)優(yōu)化的轉化,首先給出進化策略進行問題求解的框圖(圖1)和相應的步驟。

圖1 進化策略的算法流程Fig.1 Flowchart of evolutionary strategy algorithm

步驟1:產(chǎn)生初始群體。采用實數(shù)編碼,隨機產(chǎn)生初始群體,設群體規(guī)模為μ。

步驟2:計算適應度值。對所有個體將其目標函數(shù)作為適應度,適應度函數(shù)為f。當最小的適應度值f小于等于一個很接近0的值η時,算法終止。

步驟3:對群體進行排序。所有個體按照適應度值的大小進行排序,適應度值較小的個體排在適應度值較大的個體的前面。

步驟4:選擇操作。從步驟3排序完成后的群體中,選取前μ個個體。

步驟5:進化操作。將步驟4中選取的μ個個體作為父代,進行正態(tài)分布的變異遺傳操作,生成λ個子代;轉到步驟2。

用進化策略求解參數(shù)反演問題時,根據(jù)問題中參數(shù)K或B容許區(qū)間對參數(shù)采取實數(shù)編碼策略。(1)模式參數(shù)優(yōu)化。通常預報模式已知,但其參數(shù)未知,因此需要對其參數(shù)優(yōu)化。設預報模式M中有p個待參數(shù)為K={k1,k2,…,kp},反演問題轉化為一組最優(yōu)估計系數(shù)Kop,使得目標函數(shù)預報誤差f取得最小值。(2)初始場優(yōu)化。在求解初始場反問題時,設待優(yōu)化的n維初始場B={b1,b2,…,bn}。初始場反問題轉化為一組n個初始場中各變量的最優(yōu)估計,得到初始場B的最優(yōu)分析場Ba,使得目標函數(shù)預報誤差f取得最小值。

3 實驗與討論

為了進一步驗證基于進化策略的氣象學反問題求解算法有效性及其性能,進行了兩類數(shù)值模擬試驗。試驗1:是預報模式中的參數(shù)識別反演問題,采用預報模式已知、模式參數(shù)未知的線性擴散方程進行數(shù)值試驗[4];試驗2:是初始場優(yōu)化反演問題,用Lorenz-96模式進行數(shù)值試驗[12]。

3.1 一維擴散方程中的參數(shù)反演

考慮如下的問題:

取s=0.05cos(2πx-0.1t),正問題是給出擴散參數(shù)k(x)后,用差分方法求出離散解ui,m=u(iΔx,mΔt)。反問題是給出u(x,t)的觀測值ui,m,i=0,1,2,…,20,m=50,要求確定擴散參數(shù)k(x)。試驗中取一般中央差格式,Δx=0.05,Δt=0.01。

試驗采用的k(x)的真值如圖2所示,有強的間斷性。依據(jù)文獻[4]中對k(x)的估計函數(shù)形式

這樣問題轉化為對a1、a2、a3、b2、b3等5個參數(shù)的反演。

為了便于進化計算求解,構造如下目標函數(shù)

式中:uob為觀測值。此時,就可將參數(shù)反演問題轉化為下列函數(shù)優(yōu)化問題。若求得解(a1,a2,a3,b2,b3),將其帶入式(5)中,得到k＊(x),通過中央差格式求得,使得f=0(或者f為接近0的數(shù)),那么k＊(x)為擴散參數(shù)的最優(yōu)估計。

算法在MATLAB7.1上實現(xiàn),運行參數(shù)設置如下:子代規(guī)模μ=30、λ=200、pf=0.45,應用ES-(μ,λ)進化策略。為了更好地與真值kT(x)比較,評價指標采用:

統(tǒng)計參數(shù)均方誤差

經(jīng)過100代進化得到的參數(shù)k(x)的估計值:k＊(x)=0.018 500 5+0.001 285 5cos(2πx)-0.000 004 4cos(3πx)-0.000 006 1sin(2πx)+0.000 561 1sin(2πx)。相應的結果如圖2所示,實線表示待定參數(shù)k(x)的真值,帶點虛線表示通過進化策略算法求得參數(shù)k(x)的最優(yōu)估計值,此時模式預報的均方誤差為1.996×10-3。文獻[13]中用最優(yōu)k(x)估計值:k＊(x)=0.018 510 0+0.001 418 4cos(2πx)+0.000 018 8cos(3πx)-0.000 016 7sin(2πx)+0.000 384 1sin(2πx),重新計算得模式預報均方差為2.075×10-3。因此本文通過進化策略求得模式k(x)參數(shù)估計值的預報均方差更小。

圖2 一維擴散方程參數(shù)k(x)真值與反演結果Fig.2 Real value and inversion result of parameter k(x)in one-dimensional diffusion equation

3.2 Lorenz-96模式的初始場反演

Lorenz-96模式是LOREN Z和EMANU EL于1996年從動力模式中推導并簡化出的新的動力模型,已被廣泛用來驗證各類同化算法。該模型仍然具有Lorenz模型的特點,對初值極端敏感,且有很強的非線性。該模式是一個可變尺度的低階動力學系統(tǒng),模式中包含n個狀態(tài)變量X1,X2,…,Xn,均勻分布在同一緯圈上。模式的控制方程為

式中:i為循環(huán)標號;參數(shù)F=8;n一般取40。正問題是確定一組初始狀態(tài)變量X1,X2,…,Xn,通過差分格式求解T時刻的狀態(tài)變量。反問題是給出T時刻的狀態(tài)變量,要求確定初始狀態(tài)變量X1,X2,…,Xn。使用4階Runge-Kutta差分格式來求解該模式,時間步長t=0.05[13],T=mΔt,m取20。

顯然若求得解Y=(Y1,Y2,…,Yn),將其代入式(7),通過4階Runge-Kutta差分格式求出Y＊=,并且Y＊滿足f＊=0(或者f為接近0的數(shù))。那么解Y=(Y1,Y2,…,Yn)為所求初始狀態(tài)變量的最優(yōu)估計。

理想試驗方法如下:首先在預先給定理想初始場Bm,借助正問題進行預報,獲得T時刻的理想預報場;然后隨機產(chǎn)生初始場B,通過預報獲得T時刻的預報場(Y1,Y2,…,Yn);最后與理想個例的的誤差來調整和優(yōu)化初始場,進而得到最優(yōu)的初始分析場Ba和T時刻相應的預報場

算法在MA TLAB7.1上實現(xiàn),運行參數(shù)如下:子代規(guī)模μ=30、λ=200、pf=0.45,應用進化策略。作以下兩組理想試驗:(a)進化代數(shù)為4 000,10次獨立試驗,每次試驗的初始場B不同;(b)最大進化代數(shù)分別設為1 000,2 000,3 000,4 000;在不同最大代數(shù)下,分別獨立進行10次試驗。為了更好地與理想個例比較,評價指標采用:

分析場與理想初始場間均方差

T時刻的分析預報場與理想預報場間均方差

試驗a的結果如表1所示。算法獨立運行10次,每次的分析場B不同。從表1中優(yōu)化結果可以看出,模式預報均方差已經(jīng)較小,同時理想初始場Bm與分析場Ba的均方差也較小。即在不同的初始場下,算法的求解結果都非常理想,說明進化策略對初始場不敏感,算法的收斂性不受初始條件的取值范圍影響。其次,分析場均方差均大于模式預報均方差,即在分析場均方差較大的情況下(分析場均方差均大于模式預報均方差),模式預報的均方差反而較小,從而驗證了Lorenz-96模式具有非線性特征。

試驗b中,所有試驗的初始場B均相同,在不同進化代數(shù)下,算法運行10次得到的試驗結果如表2所示。從表2中試驗結果可以看出,在不同的進化代數(shù)下,經(jīng)過優(yōu)化,預報均方差的平均值、初始群體均方差的平均值都很小,再次說明此優(yōu)化算法的有效性。另外,從預報均方差的方差、分析場均方差的方差均很小看出算法非常穩(wěn)定;同時隨著進化代數(shù)的增加,預報均方差的平均值及分析場均方差的平均值越小,即算法優(yōu)化結果越好,即分析場Ba越接近理想初始場Bm。

表1 試驗a中的T時刻預報、分析場與理想個例間的均方差Table1 Mean square error of the forecast,analysis field and ideal case at time Tin experiment a

為了更好地展示預報誤差調整初始場B的過程,以最大進化代數(shù)2 000中的某次試驗為例,給出進化代數(shù)和預報誤差的關系,如圖3所示。隨著進化代數(shù)的增加,預報誤差f＊的總體趨勢逐漸減小,雖然在某些代數(shù)上存在較小的波動,但這符合進化策略采用的選擇策略特性。在算法結束時,預報誤差f＊已經(jīng)調整的較小,即進化策略反演出Lorenz-96模式的初始狀態(tài)與實際的初始狀態(tài)很接近,表明應用進化策略求解Lorenz-96模式的初始場反問題是可行的。但沒法使得預報誤差f＊為0,這恰恰也說明Lorenz-96模式存在很強的非線性。

表2 不同進化代數(shù)下的T時刻預報、分析場與理想個例間的均方差Table2 Mean square error of forecast,analysis field and ideal case at time Tin different evolution generations

圖3 預報誤差f與進化代數(shù)之間的關系Fig.3 Relationship be tween forecast errorf and evolution generations

4 結論與討論

在研究氣象學上的天氣數(shù)值預報時,提出了一種基于進化策略的氣象學反問題求解算法,并借助預報模式中的參數(shù)識別和初始場反問題,通過問題轉化變成等價函數(shù)優(yōu)化問題,進而采用進化策略求解。為驗證方法的有效性,通過一維擴散方程和Lorenz-96預報模式中的反問題進行理想數(shù)值試驗,實驗結果表明應用進化策略求解預報模式中的參數(shù)識別問題和初始場反問題是可行的,算法都能求解到理想的結果。但對于初始場中的變量協(xié)調問題,可轉化為帶約束優(yōu)化問題,是下一步的研究方向。

當然,依據(jù)現(xiàn)有的計算機發(fā)展水平,暫且無法將進化策略直接用于當前業(yè)務預報模式中(如W RF、MM5等)。主要是因為:一是進化算法的優(yōu)化精髓是“生成—測試”方法,即需進行頻繁評估;二是由于業(yè)務模式非常復雜,在其上進行一次評估是非常耗時的,即無法滿足業(yè)務對時間的要求。我們也堅信隨著計算機硬件和并行算法的不斷發(fā)展,在不遠的將來進化算法也將為推動氣象研究與發(fā)展作出努力。

[1] 任宏利,丑紀范.數(shù)值模式的預報策略和方法研究進展[J].地球科學進展,2007,22(4):376-385.

[2] Kalnay E.A tmospheric modeling,data assim ilation and predictability[M].Cambridge:Cambridge University Press,2003.

[3] 郜吉東,丑紀范.數(shù)值天氣預報中兩類反問題及一種數(shù)值解法[J].氣象學報,1994,52(2):129-137.

[4] 邱崇踐,丑紀范.預報模式的參數(shù)優(yōu)化方法[J].中國科學,1990,2(2):218-224.

[5] 吳志健.演化優(yōu)化及其在微分方程反問題的應用[D].武漢:武漢大學計算機學院,2004.

[6] Coello C A C.Twenty years of evolutionary multi-objective optimization:A historical view of the field[J].IEEE Computational Intelligence Magazine,2006,1(1):28-36.

[7] Thomas P R,Xin Yao.Search biases in constrained evolutionary optimization[J].IEEE Transactions on Systems,2005,5(2):233-243.

[8] 常新功,李敏強,寇紀淞.基于進化算法的圖形數(shù)據(jù)模式發(fā)現(xiàn)[J].模式識別與人工智能,2008,21(1):116-121.

[9] 馬清亮,胡昌華.多目標進化算法及其在控制領域中的應用綜述[J].控制與決策,2006,21(5):481-486.

[10] 李守巨,劉迎曦,孫偉.智能計算與參數(shù)反演[M].北京:科學出版社,2008.

[11] 韓靜.利用星載超光譜儀器觀測資料進行大氣參數(shù)反演的試驗[D].南京:南京信息工程大學遙感學院,2008.

[12] 王文強.基于Lorenz-96模式的集合平方根濾波同化方法研究[D].蘭州:蘭州大學大氣科學學院,2007.

[13] 熊盛武,李元香,康立山,等.拋物型方程的演化參數(shù)識別方法[J].計算物理,2000,17(5):511-517.