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

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

0 引言

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

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

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

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

1 問題描述

1.1 微分方程反問題

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

1.2 微分方程反問題的分類

設(shè)研究的微分方程一般形式如下:

微分方程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分別為微分算子、初始算子、邊界算子和附加算子。當(dāng)上述這些已知量中有一個變?yōu)槲粗獣r,就構(gòu)成了偏微分方程的反問題。

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

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

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

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

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

2 應(yīng)用進(jìn)化策略求解微分方程反問題

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

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

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

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

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

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

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

2.2 進(jìn)化策略

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

圖1 進(jìn)化策略的算法流程Fig.1 Flowchart of evolutionary strategy algorithm

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

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

步驟3:對群體進(jìn)行排序。所有個體按照適應(yīng)度值的大小進(jìn)行排序,適應(yīng)度值較小的個體排在適應(yīng)度值較大的個體的前面。

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

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

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

3 實驗與討論

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

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

考慮如下的問題:

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

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

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

為了便于進(jìn)化計算求解,構(gòu)造如下目標(biāo)函數(shù)

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

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

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

經(jīng)過100代進(jìn)化得到的參數(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)。相應(yīng)的結(jié)果如圖2所示,實線表示待定參數(shù)k(x)的真值,帶點(diǎn)虛線表示通過進(jìn)化策略算法求得參數(shù)k(x)的最優(yōu)估計值,此時模式預(yù)報的均方誤差為1.996×10-3。文獻(xiàn)[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),重新計算得模式預(yù)報均方差為2.075×10-3。因此本文通過進(jìn)化策略求得模式k(x)參數(shù)估計值的預(yù)報均方差更小。

圖2 一維擴(kuò)散方程參數(shù)k(x)真值與反演結(jié)果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年從動力模式中推導(dǎo)并簡化出的新的動力模型,已被廣泛用來驗證各類同化算法。該模型仍然具有Lorenz模型的特點(diǎn),對初值極端敏感,且有很強(qiáng)的非線性。該模式是一個可變尺度的低階動力學(xué)系統(tǒng),模式中包含n個狀態(tài)變量X1,X2,…,Xn,均勻分布在同一緯圈上。模式的控制方程為

式中:i為循環(huán)標(biāo)號;參數(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)估計。

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

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

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

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

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

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

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

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

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

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

4 結(jié)論與討論

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

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

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

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

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

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

[5] 吳志健.演化優(yōu)化及其在微分方程反問題的應(yīng)用[D].武漢:武漢大學(xué)計算機(jī)學(xué)院,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] 常新功,李敏強(qiáng),寇紀(jì)淞.基于進(jìn)化算法的圖形數(shù)據(jù)模式發(fā)現(xiàn)[J].模式識別與人工智能,2008,21(1):116-121.

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

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

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

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

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