顏冰 黃思訓(xùn)2)? 馮徑

1)(國(guó)防科技大學(xué)氣象海洋學(xué)院,南京 211101)

2)(上海財(cái)經(jīng)大學(xué),計(jì)算科學(xué)與金融數(shù)據(jù)研究中心,上海 200433)

(2018年5月24日收到；2018年7月10日收到修改稿)

1 引 言

行星大氣邊界層的研究對(duì)于氣候模式和數(shù)值預(yù)報(bào)模式的發(fā)展很重要,邊界層模式的選擇非常關(guān)鍵[1].完整的邊界層模式過于復(fù)雜不利于計(jì)算與分析,簡(jiǎn)化的Ekman模式過于簡(jiǎn)單無法描述邊界層的非線性特征.Tan和Wang[2]提出的廣義Ekman動(dòng)量近似模式既保留了慣性項(xiàng)的作用,又考慮了湍流黏性系數(shù)隨高度的變化,是一個(gè)研究大氣邊界層實(shí)用且有效的模式.在該模式中,一個(gè)非常重要的參數(shù)就是湍流黏性系數(shù),它是反映湍流效應(yīng)的參數(shù)[3].通常,直接測(cè)量出準(zhǔn)確的湍流黏性系數(shù)較為困難,所以人們希望得到不確定性較小的湍流黏性系數(shù).

根據(jù)湍流黏性系數(shù)的分布計(jì)算出風(fēng)速的分布,叫做“正向不確定量化”或者 “預(yù)報(bào)模式”[4].實(shí)際上湍流黏性系數(shù)是確定的,但是由于觀測(cè)誤差等原因無法準(zhǔn)確獲知其準(zhǔn)確值,僅能得到大致的分布,所以必須收集更多的信息來縮小分布的范圍.通過風(fēng)速的觀測(cè)反演降低湍流黏性系數(shù)的不確定性,叫做“逆向不確定量化”或者 “逆模式”[5].不確定性問題的研究比確定性問題困難,同時(shí),因?yàn)榉磫栴}是不適定問題,這意味著解可能不存在、解存在但是多解或解不穩(wěn)定.反問題的研究比正問題的研究困難很多[6],所以逆向不確定問題的研究更加困難,而且由于觀測(cè)數(shù)據(jù)存在著不可避免的誤差,觀測(cè)中包含的信息通常不足以確定所有的不確定參數(shù),所以參數(shù)的不確定性不能被消除,只能降低.

反問題一般根據(jù)貝葉斯理論來求解,利用給定的參數(shù)先驗(yàn)分布結(jié)合模式輸出量計(jì)算得到參數(shù)的后驗(yàn)分布.但是它在大部分情況下都不能被解析地表達(dá),所以必須通過采樣的方法來表達(dá)后驗(yàn)分布的不確定性,常用的方法有馬爾可夫鏈蒙特卡羅抽樣方法(Markov chain Monte Carlo,MCMC),基于Metropolis-Hastings算法[7],Gibbs采樣器[8],粒子濾波[9].但是這些采樣的方法需要大量的樣本,每一個(gè)樣本均需要進(jìn)行一次模式的計(jì)算,所以占用了大量的計(jì)算資源.

一種常用的能夠解析計(jì)算后驗(yàn)分布的方法叫做Kalman濾波[10],它的要求是模式必須是線性的,而且先驗(yàn)分布與觀測(cè)誤差分布必須是正態(tài)分布.這種情況下模式參數(shù)的后驗(yàn)分布同樣是正態(tài)分布,它的特征能夠用前兩階矩即均值與方差來表征.但是實(shí)際大部分大氣海洋數(shù)值模型是非線性的,而且觀測(cè)模型也是非線性的.對(duì)于非線性系統(tǒng),可以用Kalman濾波得到一個(gè)較為準(zhǔn)確的近似解,其難點(diǎn)在于在耗費(fèi)較少計(jì)算資源的情況下怎樣由參數(shù)的分布計(jì)算得到模式輸出的分布.擴(kuò)展Kalman濾波方法(extended Kalman filter,EKF)是Kalman濾波的近似,將非線性模式在分析值處進(jìn)行Taylor展開到一階,略去高階項(xiàng),再進(jìn)行Kalman濾波過程[11].但是在計(jì)算過程中由于非線性會(huì)導(dǎo)致協(xié)方差矩陣的不穩(wěn)定性,而且計(jì)算量大.集合Kalman濾波方法(ensemble Kalman filter,EnKF)首先利用蒙特卡羅法來求解不確定性在非線性系統(tǒng)中的傳播[12,13],然后計(jì)算得到Kalman增益矩陣.Evensen[13]最早把EnKF用到數(shù)據(jù)同化中,他利用概率論上一些理論與簡(jiǎn)單數(shù)值實(shí)驗(yàn)說明該方法的可行性.由于EnKF便于實(shí)施,故在很多領(lǐng)域都得到了廣泛的應(yīng)用.但是由于EnKF基于蒙特卡羅方法,所以其準(zhǔn)確性就與樣本的大小密切相關(guān),而大模式樣本需要非常多的計(jì)算資源.Kalman濾波類的方法通常是采用遞歸的方式來進(jìn)行,對(duì)于每個(gè)數(shù)據(jù)序列上的數(shù)據(jù)均進(jìn)行一次濾波計(jì)算,這樣前一次的后驗(yàn)分布作為后一次的先驗(yàn)分布,這種過程被稱為“順序數(shù)據(jù)同化”.

混沌多項(xiàng)式展開(polynomial chaos expansion,PCE)可以作為蒙特卡羅法的代替,其本質(zhì)就是利用一組正交隨機(jī)多項(xiàng)式展開來表達(dá)不確定性[14,15].根據(jù)Askey方案[16],當(dāng)分布已知的情況下,最優(yōu)的正交基函數(shù)就是確定的,所需要求解的就是待定的多項(xiàng)式系數(shù).當(dāng)系數(shù)求解出來后,就會(huì)很快求出隨機(jī)變量的統(tǒng)計(jì)矩[17].所以混沌多項(xiàng)式展開方法的關(guān)鍵就是求解多項(xiàng)式系數(shù).有很多方法可以求解此系數(shù),例如隨機(jī)Galerkin方法[15],回歸法[18],概率配點(diǎn)法[19].本文中利用回歸法來求解系數(shù),它是一種非浸入式方法,不需要對(duì)原模式計(jì)算代碼進(jìn)行更改,只需要進(jìn)行少量的采樣,然后求解方程組即可求得系數(shù).該方法易于理解,便于實(shí)施,而且在消耗較少資源的情況下能達(dá)到較高的精度.模式輸入與輸出變量采用混沌多項(xiàng)式表達(dá),利用基函數(shù)的正交性求出輸出變量的統(tǒng)計(jì)矩,避免了使用蒙特卡羅法,能夠節(jié)省計(jì)算資源.這種方法就被稱為基于混沌多項(xiàng)式的集合Kalman濾波方法(polynomial chaos-ensemble Kalman filter,PC-EnKF)[20].

本文試圖利用PC-EnKF求解觀測(cè)值為風(fēng)速的情況下湍流黏性系數(shù)的后驗(yàn)分布,以降低系數(shù)不確定性和提供更可靠的范圍.而且,還可以根據(jù)先驗(yàn)的風(fēng)速分布確定出有效觀測(cè)的區(qū)域,為觀測(cè)點(diǎn)位置的選擇提供強(qiáng)有力的指導(dǎo).由于實(shí)際觀測(cè)數(shù)據(jù)難以獲得,所以本文采用數(shù)值實(shí)驗(yàn)的方式驗(yàn)證該方法的有效性和可靠性.本文第2部分介紹隨機(jī)廣義Ekman動(dòng)量近似模式以及預(yù)報(bào)模式的求解方法；第3部分介紹逆模式以及Kalman濾波類算法；第4部分介紹數(shù)值實(shí)驗(yàn)及其分析；第5部分給出總結(jié)與討論.

2 隨機(jī)廣義Ekman動(dòng)量近似模式

2.1 模 式

本文采用Tan和Wang[2]提出的廣義Ekman動(dòng)量近似模式(general Ekman momentum approximation model,GEM),該模式簡(jiǎn)單且同時(shí)保留了系統(tǒng)的非線性特征.將模式中湍流黏性系數(shù)由確定性的參數(shù)變?yōu)榉夏撤N分布的不確定性參數(shù),此時(shí)該模式演變成廣義Ekman動(dòng)量近似的隨機(jī)邊界層模式(stochastic general Ekman momentum approximation model,SGEM).模型如下:

其中湍流黏性系數(shù)K是隨高度z變化的隨機(jī)變量,ξ代表的是隨機(jī)變量,(u,v)代表(x,y)方向上的風(fēng)速分量；

(2)式中,f代表Coriolis參數(shù),在本文中假定為常數(shù)；(ue,ve)為經(jīng)典Ekman模式的解,(ug,vg)為地轉(zhuǎn)風(fēng).方程(1)的邊界條件為:u|z=0=0,v|z=0=0；u|z=Hb=uT,v|z=Hb=vT.其中Hb代表邊界層的厚度,在本文中假設(shè)Hb=1500 m；(uT,vT)代表邊界層頂部的風(fēng)速,可由地轉(zhuǎn)動(dòng)量近似模型求得.本文在空間、時(shí)間方向上將方程(1)離散化,然后直接進(jìn)行數(shù)值計(jì)算.

2.2 SGEM正模型求解

對(duì)于此隨機(jī)模式的解法,本文擬采用混沌多項(xiàng)式展開的方法.根據(jù)Askey方案[16],若湍流黏性系數(shù)K的分布已知,其對(duì)應(yīng)的最優(yōu)混沌多項(xiàng)式基函數(shù)Ψ也就確定下來了.然后將輸入隨機(jī)變量K與輸出隨機(jī)量(u,v)在此基函數(shù)上展開,不同的是K的系數(shù)是已知的,而(u,v)的系數(shù)是需要求解的.所以問題由求解(u,v)的分布變?yōu)榍蠼馄涠囗?xiàng)式展開的系數(shù).本文中采用回歸法求解展開系數(shù),此方法簡(jiǎn)單有效,易于理解,在使用較少計(jì)算資源的情況下也能有不錯(cuò)的效果.主要的思想就是對(duì)隨機(jī)變量采樣(樣本數(shù)應(yīng)該大于兩倍的混沌多項(xiàng)式的項(xiàng)數(shù))[17],然后計(jì)算得到解的樣本集和多項(xiàng)式基函數(shù)的樣本集,最后求解超定方程組,結(jié)果為多項(xiàng)式展開的系數(shù).將SGEM模式記為U=g(K),輸入隨機(jī)變量K為湍流黏性系數(shù),輸出隨機(jī)變量TΨi(ξ)為多項(xiàng)式基函數(shù),ξ =[ξ1,ξ2,...,ξd]；P+1為多項(xiàng)式的個(gè)數(shù),它與隨機(jī)變量的維數(shù)d和截?cái)嗟碾A數(shù)r有關(guān),即P+1=(d+r)!/(d!r!)[4].求解流程如圖1所示.

3 逆模式及Kalman濾波類算法

若模式參數(shù)K是未知的,但我們可以根據(jù)專家猜值或者一些直接的但精度有限的觀測(cè)得知其先驗(yàn)的概率密度函數(shù)p(K),風(fēng)速U可以通過觀測(cè)得到.為了降低參數(shù)的不確定性,得到一個(gè)更加接近真值的估計(jì)值,可以利用帶誤差的觀測(cè)U?.U?=U+eU,誤差eU符合另一概率密度函數(shù).求解逆模式的目的就是估計(jì)出參數(shù)K在已知觀測(cè)U?的情況下的后驗(yàn)概率密度p(K|U?).

圖1 SGEM預(yù)報(bào)模式求解流程圖Fig.1.the flow chart for solving the SGEM.

3.1 Kalman濾波

Kalman濾波能給出準(zhǔn)確的后驗(yàn)概率密度函數(shù),但它有著較高的限制條件:1)模式必須是線性的；2)參數(shù)的先驗(yàn)密度函數(shù)p(K)與誤差eU均服從正態(tài)分布.此時(shí),參數(shù)的后驗(yàn)概率密度函數(shù)同樣服從正態(tài)分布,其特征可以用均值μK|U?與協(xié)方差CKK|U?來表達(dá).其求解過程如下.

第一步 預(yù)測(cè),利用參數(shù)K先驗(yàn)分布通過線性隨機(jī)模式求出輸出量U的分布,又叫做預(yù)報(bào)模式求解過程.由于模式是線性的且均服從正態(tài)分布,所以很容易得到μU,CUU,CUK.μU代表U的先驗(yàn)/預(yù)測(cè)均值,CUK是U和K的先驗(yàn)協(xié)方差,CUU是向量U的先驗(yàn)/預(yù)測(cè)協(xié)方差.然后計(jì)算得到Kalman增益矩陣G.G=CKU(CUU+R)?1,式中,CKU=(CUK)T,R是誤差eU的協(xié)方差矩陣,通常是預(yù)先給定的.

第二步 分析.通過觀測(cè)值來調(diào)整參數(shù)K的分布,求出后驗(yàn)分布,具體來說就是求出后驗(yàn)均值與后驗(yàn)方差,可以通過下列公式得到:

其中μK和CKK分別代表參數(shù)K的先驗(yàn)的均值和協(xié)方差.

通常,觀測(cè)數(shù)據(jù)是一個(gè)序列,所以Kalman濾波也會(huì)通過遞歸進(jìn)行,每一個(gè)觀測(cè)數(shù)據(jù)均會(huì)進(jìn)行一次Kalman濾波過程.也就是說,根據(jù)前一次觀測(cè)數(shù)據(jù)而求得后驗(yàn)分布會(huì)作為后一次Kalman濾波的先驗(yàn)分布.

3.2 EnKF

Kalman濾波方法只適用于線性模式,而大氣海洋科學(xué)中大部分模式都是非線性的.對(duì)于非線性不太強(qiáng)的模式,我們可以用EnKF方法給出一個(gè)近似解,得出的后驗(yàn)分布不一定符合正態(tài)分布,但也是單峰型的,可以用前兩階矩(均值和方差)來表示其特征.對(duì)于非線性模式,通常采用蒙特卡羅法來求解輸出變量的先驗(yàn)/預(yù)測(cè)分布,這就是EnKF.首先根據(jù)輸入?yún)?shù)的先驗(yàn)分布采樣,得到一組參數(shù)的樣本集,然后對(duì)每一個(gè)樣本進(jìn)行確定預(yù)報(bào)模式的計(jì)算,得到輸出變量的樣本集.

在輸出變量樣本集的基礎(chǔ)上可以計(jì)算得到均值、協(xié)方差以及與輸入?yún)?shù)的協(xié)方差.

然后是分析步,由于是采樣的方法,需要對(duì)每一個(gè)樣本進(jìn)行更新:

之后同樣進(jìn)行遞歸的過程.可以看到在(6)式中需要對(duì)觀測(cè)U?根據(jù)給定的協(xié)方差進(jìn)行擾動(dòng)后采樣得到Ui?,這就會(huì)產(chǎn)生不必要的采樣誤差.Whitaker和 Hamill[21]提出了集合均方根濾波方法(ensemble square root filter,EnSRF),此方法可以使用確定性的觀測(cè)數(shù)據(jù),這樣就可以避免在觀測(cè)時(shí)產(chǎn)生采樣的誤差.這種方法將參數(shù)K的均值和擾動(dòng)分別進(jìn)行更新.均值的更新使用標(biāo)準(zhǔn)的Kalman增益矩陣,而擾動(dòng)的更新使用修正的增益矩陣:

在本文中,我們EnSRF方法來進(jìn)行反演的過程.

3.3 PC-EnKF

EnKF使用蒙特卡羅法進(jìn)行輸出變量的不確定量化,而蒙特卡羅法的收斂速度慢且需要足夠大的樣本來保證結(jié)果的準(zhǔn)確性,這就造成了巨大的計(jì)算負(fù)擔(dān),特別是對(duì)于大型非線性模式.為了在不顯著增加計(jì)算量的情況下增加逆模式的求解精度,本文采用EnKF的另外一個(gè)變種:基于混沌多項(xiàng)式的集合Kalman濾波方法(PC-EnKF).此方法類似于EnKF,只是在預(yù)報(bào)模式中采用混沌多項(xiàng)式展開來代替蒙特卡羅法進(jìn)行不確定性量化.由于混沌多項(xiàng)式展開的關(guān)鍵在于多項(xiàng)式的系數(shù),所以將EnKF中樣本的更新改變?yōu)槎囗?xiàng)式系數(shù)的更新.具體實(shí)施步驟如下.

第一步 預(yù)測(cè). 首先根據(jù)參數(shù)K的先驗(yàn)分布找出其最優(yōu)的多項(xiàng)式基函數(shù)Ψi(ξ),ξ=[ξ1,ξ2,...,ξd]T是由一組給定的概率密度函數(shù)的相互獨(dú)立的隨機(jī)變量組成的,例如標(biāo)準(zhǔn)正態(tài)分布、均勻分布等.然后在此基函數(shù)上展開得到:由于其先驗(yàn)分布是已知的,所以系數(shù)ki也是已知的.然后利用上文說到的預(yù)報(bào)模式求解方法,求出中Ui是多項(xiàng)式展開系數(shù).在本文中,假設(shè)參數(shù)服從正態(tài)分布,其對(duì)應(yīng)的最優(yōu)多項(xiàng)式為Hermite多項(xiàng)式[4,22].混沌多項(xiàng)式都有著標(biāo)準(zhǔn)正交的特性,E(Ψi(ξ)Ψj(ξ))= δij,δij是Kronecker delta函數(shù).

結(jié)合誤差eU的協(xié)方差R,而且注意到CKU=(CUK)T,可以計(jì)算出Kalman增益和修正的增益.

第二步 分析.將觀測(cè)數(shù)據(jù)同化進(jìn)來,由此調(diào)整得到參數(shù)K的后驗(yàn)分布的多項(xiàng)式展開系數(shù).本文中我們采用EnSRF策略進(jìn)行更新,所以分為均值更新與擾動(dòng)更新:

圖2 利用PC-EnKF減少參數(shù)不確定性求解流程圖Fig.2.the flow chart for reducing the uncertainty by PC-EnKF.

如有多次的觀測(cè),可利用遞歸的方式進(jìn)行計(jì)算,根據(jù)前一次觀測(cè)數(shù)據(jù)而求得后驗(yàn)分布會(huì)作為后一次Kalman濾波的先驗(yàn)分布.

后驗(yàn)分布的均值與方差容易求得:

其流程如圖2所示.

4 數(shù)值實(shí)驗(yàn)

為了驗(yàn)證PC-EnKF在反演SGEM中湍流黏性系數(shù)K不確定性時(shí)的有效性.本文選取一類典型的氣旋性地轉(zhuǎn)切變流進(jìn)行實(shí)驗(yàn),其表達(dá)式為[2]:

其中u0=20 m.s?1；α為常數(shù),取α =0.4f,而f=10?4s?1.由于實(shí)際觀測(cè)數(shù)據(jù)較難獲得,本文決定采用模擬觀測(cè)數(shù)據(jù)進(jìn)行數(shù)值實(shí)驗(yàn),在后續(xù)的研究中再使用實(shí)際觀測(cè)數(shù)據(jù).

我們假設(shè)反演兩種類型下的湍流黏性系數(shù)；第一種,湍流黏性系數(shù)是常數(shù),不隨高度變化；第二種,湍流黏性系數(shù)是隨高度變化的量.

何良諸打車，趕到北大坎市公安局。市局對(duì)姍姍來遲的何良諸，已經(jīng)不感興趣。礦工們上訪、鬧事，事出有因，領(lǐng)導(dǎo)一直下不了決心，執(zhí)法部門不能動(dòng)手。市局原定以倒賣稀世文物、查抄琥珀銘文為由，控制趙集、小勺等人。市局將方案報(bào)省后，既熟悉北大坎礦區(qū)，又是文物專家的何良諸被派下來。攔截國(guó)線火車的事件發(fā)生后，市領(lǐng)導(dǎo)抓住把柄，有膽量拍板了。那個(gè)晚上，警察出動(dòng)，拘捕幾個(gè)人，敲山震虎，天下太平了。

4.1 實(shí)驗(yàn)一

為了避免湍流黏性系數(shù)出現(xiàn)負(fù)值的情況,我們假設(shè)π=log(K),π服從正態(tài)分布.根據(jù)湍流黏性系數(shù)取值的大致范圍確定先驗(yàn)分布π～N(2,0.42),然后按照先驗(yàn)分布隨機(jī)采樣得到一個(gè)樣本作為參考的真值參數(shù),然后根據(jù)確定性的GEM模式計(jì)算得到相應(yīng)的風(fēng)速,在此基礎(chǔ)上疊加符合正態(tài)分布的觀測(cè)誤差R～N(0,0.22),N代表正態(tài)分布,各點(diǎn)的誤差假定是相互獨(dú)立的,由此生成模擬觀測(cè)數(shù)據(jù).觀測(cè)數(shù)據(jù)隨高度變化U?(l)=U?(zl),l=1,2,...,NH,NH為垂直方向格點(diǎn)數(shù),在每個(gè)點(diǎn)上均進(jìn)行一次Kalman濾波,這樣就形成一個(gè)遞歸的過程.我們將垂直方向分成51個(gè)格點(diǎn),這樣就會(huì)有49個(gè)Kalman濾波過程(因?yàn)橐コ齼蓚€(gè)固定邊界條件).按照?qǐng)D2所示的流程圖進(jìn)行反演,為了直觀地顯示遞歸過程中的不確定性降低,每進(jìn)行十六次循環(huán),繪制一次模式輸入量π和對(duì)應(yīng)的模式輸出量U的后驗(yàn)分布的均值和方差圖形.

圖3 模式參數(shù)π的置信區(qū)間(均值±標(biāo)準(zhǔn)差) 橫坐標(biāo)代表經(jīng)過的迭代步數(shù),0步代表并未經(jīng)過同化的先驗(yàn)分布,其余為經(jīng)過同化的后驗(yàn)分布；縱坐標(biāo)代表π的取值；星號(hào)代表π的參考真值,方塊代表π的后驗(yàn)分布的均值；注意:在(a)中,由于后驗(yàn)分布的不確定性很小,而(a)的縱坐標(biāo)范圍過大不能顯示其中的細(xì)節(jié),故將后驗(yàn)分布的置信區(qū)間放大放在(b)中Fig.3.The estimated confidence interval of model parameter π (mean ± standard deviation).The horizontal coordinate is the iterative number,where n=0 represents the prior distribution without assimilation,and the others are the posterior distribution after assimilation.The vertical coordinate is the value of π.The asterisk represents the reference value,the square represents the mean of posterior distribution.Note that the uncertainty of posterior distribution is small,while the range of vertical coordinate of Fig.(a)is too large to display details of posterior distribution.So,these details are shown in Fig.(b).

圖4 風(fēng)速Ekman螺線的觀測(cè)數(shù)據(jù)以及風(fēng)速分布示意圖 灰色區(qū)域是分布的置信區(qū)間(均值±標(biāo)準(zhǔn)差),虛線代表Ekman螺線分布的均值,空心圓點(diǎn)代表觀測(cè)點(diǎn),也是數(shù)據(jù)同化點(diǎn)Fig.4.The wind observations of Ekman spirals and the distribution of wind.The grey shaded area is the confidence interval of wind speed(mean±standard deviation),the dotted line is the mean of wind distribution,the circular point is the observation point and also the data assimilation point.

從圖3(a)可以看出,PC-EnKF能有效利用風(fēng)速的觀測(cè)信息進(jìn)行反演,顯著降低參數(shù)的不確定性.從圖3(b)可以看出,隨著同化過程的進(jìn)行,參數(shù)不確定性在逐步減小且其后驗(yàn)分布的均值越來越靠近參考真值.

圖4(a)中Ekman螺線先驗(yàn)分布的不確定性較大,在經(jīng)過同化步后(圖4(bd)),風(fēng)速的不確定性明顯降低,而且與觀測(cè)數(shù)據(jù)的擬合程度較高,說明了PC-EnKF在降低參數(shù)不確定性方面是有效的.經(jīng)過16個(gè)同化步之后預(yù)測(cè)的Ekman螺線與觀測(cè)已經(jīng)比較接近,也對(duì)應(yīng)著圖3(b)中經(jīng)過16步的同化過程后模式參數(shù)的后驗(yàn)分布小的不確定性.

本文在風(fēng)速先驗(yàn)分布中,不確定性較大的區(qū)域與不確定性小的區(qū)域各選取了10個(gè)觀測(cè)點(diǎn),在此基礎(chǔ)上分別進(jìn)行反演,比較其反演結(jié)果.

根據(jù)圖5所示的觀測(cè)點(diǎn)上的觀測(cè)數(shù)據(jù)分別進(jìn)行反演,可以分別求得兩種觀測(cè)下的參數(shù)后驗(yàn)分布.

從圖6中可以明顯地看出,相比于b區(qū)域,a區(qū)域觀測(cè)的反演效果更好,均值更靠近參考真值且標(biāo)準(zhǔn)差更小.因?yàn)榉囱輰儆诜磫栴},反問題的一個(gè)重要特點(diǎn)就是不穩(wěn)定,不穩(wěn)定即模式輸出值小的變化就會(huì)引起參數(shù)的反演結(jié)果很大的誤差,對(duì)應(yīng)到不確定性分析中就是模式參數(shù)的先驗(yàn)分布不確定性較大,而對(duì)應(yīng)的模式輸出的先驗(yàn)分布不確定性較小.由一個(gè)小范圍的分布去反演一個(gè)大范圍的分布就會(huì)引起反演結(jié)果的不穩(wěn)定.由此我們認(rèn)為在模式參數(shù)先驗(yàn)分布不變的情況下,選取輸出值先驗(yàn)分布不確定性較大的區(qū)域的觀測(cè)進(jìn)行反演會(huì)更加有效率,圖6的結(jié)果也證明了我們的觀點(diǎn).這個(gè)發(fā)現(xiàn)對(duì)于風(fēng)速觀測(cè)點(diǎn)的位置的選擇具有指導(dǎo)意義.

圖5 選取的觀測(cè)點(diǎn)示意圖 灰色區(qū)域是風(fēng)速先驗(yàn)分布的置信區(qū)間(均值±標(biāo)準(zhǔn)差),星號(hào)代表觀測(cè)區(qū)域a(不確定性大的區(qū)域)的觀測(cè),圓圈代表觀測(cè)區(qū)域b(不確定性小的區(qū)域)的觀測(cè)Fig.5.The selected observation points.The grey shaded area is the confidence interval of wind speed(mean±standard deviation).The star represents observation in area a(heavy uncertainty),the circle represents observation in area b(small uncertainty).

圖6 由不同區(qū)域的觀測(cè)反演的結(jié)果 橫坐標(biāo)代表不同的區(qū)域,縱坐標(biāo)的模式參數(shù)π的取值,星號(hào)代表π的參考真值,方塊代表π的后驗(yàn)分布的均值Fig.6.The result retrieved by different observations.The horizontal coordinate represents the area,the vertical coordinate is the value of π.The asterisk represents the reference,the square represents the mean of posterior distribution.

4.2 實(shí)驗(yàn)二

考慮湍流黏性系數(shù)一般是隨高度變化的量,我們假設(shè)每個(gè)點(diǎn)的湍流黏性系數(shù)都是一個(gè)隨機(jī)變量.因此系統(tǒng)成為高維的,混沌多項(xiàng)式的項(xiàng)數(shù)會(huì)急劇增加,計(jì)算量也會(huì)急劇增加.在本文中,我們假設(shè)垂直方向上格點(diǎn)的湍流黏性系數(shù)是存在相關(guān)性.K(z,ξ)是模式參數(shù),假設(shè)其先驗(yàn)分布為對(duì)數(shù)正態(tài)分布,假設(shè)A=log(K),A符合正態(tài)分布,E(A)=2.08,Cov(A(z1),A(z2))=exp(?|z1?z2|/c),a=25,c=0.004. 我們?cè)诖嘶A(chǔ)上隨機(jī)生成參考真值,然后經(jīng)過確定性模式得到模式輸出值,再加上符合觀測(cè)誤差協(xié)方差R～N(0,0.22)的誤差變成了帶誤差的觀測(cè)值,N代表正態(tài)分布,且假設(shè)各點(diǎn)的觀測(cè)誤差獨(dú)立.

在利用PC-EnKF進(jìn)行反演之前,必須將模式參數(shù)A的先驗(yàn)分布表達(dá)成混沌多項(xiàng)式展開的形式,利用K-L分解(Karhunen-Loève expansion)[4],由于A的先驗(yàn)均值與協(xié)方差函數(shù)已知,可以將其近似表達(dá)為

詳細(xì)的步驟可以參閱文獻(xiàn)[4,22].

圖7 模式參數(shù)A的分布示意 橫坐標(biāo)代表高度,縱坐標(biāo)是參數(shù)的取值；灰色區(qū)域是分布的置信區(qū)間(均值±標(biāo)準(zhǔn)差),點(diǎn)代表模式參數(shù)A的參考真值,虛線代表的是分布的均值Fig.7.The distribution of model parameter A.The horizontal coordinate represents the height,the vertical coordinate is the value of model parameter.The grey shaded area is the confidence interval of A(mean±standard deviation),the solid point is the reference value of A,the dotted line is the mean of distribution.

圖8 風(fēng)速Ekman螺線的觀測(cè)數(shù)據(jù)以及風(fēng)速分布示意圖 灰色區(qū)域是分布的置信區(qū)間(均值±標(biāo)準(zhǔn)差),虛線代表Ekman螺線分布的均值,空心圓點(diǎn)代表觀測(cè)點(diǎn),也是數(shù)據(jù)同化點(diǎn)Fig.8.The wind observations of Ekman spirals and the distribution of wind.The grey shaded area is the confidence interval of wind speed(mean±standard deviation),the dotted line is the mean of wind distribution,the circular point is the observation point and data assimilation point.

在得到先驗(yàn)分布的多項(xiàng)式展開之后,就可以實(shí)施PC-EnKF反演過程.與實(shí)驗(yàn)一類似,實(shí)驗(yàn)二同樣是每隔16次繪制一次模式輸入量π和對(duì)應(yīng)的模式輸出量U的后驗(yàn)分布的均值和方差的圖形.

從圖7(a)可以看出參數(shù)的先驗(yàn)分布不確定性是相當(dāng)大的,參考真值與分布的均值相去甚遠(yuǎn).經(jīng)過同化過程后,不確定性有了明顯的減小,特別是均值與參考真值已經(jīng)很接近了.

圖8(a)中Ekman螺線先驗(yàn)分布的不確定性不是很大.在經(jīng)過同化步后(圖8(b)—(d)),Ekman螺線的不確定性有所降低,與觀測(cè)數(shù)據(jù)的擬合程度高,說明了PC-EnKF在降低參數(shù)不確定性方面是有效的.在實(shí)驗(yàn)二中我們同樣選取了不確定性大與不確定性小兩個(gè)區(qū)域的觀測(cè),并分別進(jìn)行反演.觀測(cè)點(diǎn)的位置如圖9所示.

可以很明顯地看出,圖10(a)中的反演效果更好,具有更小的不確定性,進(jìn)一步縮小了參數(shù)分布的區(qū)間.這也佐證了我們的觀點(diǎn),在不確定性大的區(qū)域的觀測(cè)數(shù)據(jù)在反演時(shí)具有更好的效果.而且實(shí)驗(yàn)二中風(fēng)速的先驗(yàn)分布(圖9)的不確定性明顯比實(shí)驗(yàn)一(圖5)的小,根據(jù)我們上文中的分析,實(shí)驗(yàn)二中模式參數(shù)的反演的不穩(wěn)定性更強(qiáng),反演難度大于實(shí)驗(yàn)一.比較圖7與圖3,可以看出實(shí)驗(yàn)二的反演結(jié)果也較實(shí)驗(yàn)一差.

圖9 選取的觀測(cè)點(diǎn)示意圖 灰色區(qū)域是風(fēng)速先驗(yàn)分布的置信區(qū)間(均值±標(biāo)準(zhǔn)差).星號(hào)代表觀測(cè)區(qū)域a(不確定性大的區(qū)域)的觀測(cè),圓圈代表觀測(cè)區(qū)域b(不確定性小的區(qū)域)的觀測(cè)Fig.9.The selected observation points.The grey shaded area is the confidence interval of wind speed(mean±standard deviation).The star represents observation in area a(heavy uncertainty),the circle represents observation in area b(small uncertainty).

圖10 由不同區(qū)域觀測(cè)反演參數(shù)A的結(jié)果 橫坐標(biāo)代表高度,縱坐標(biāo)是參數(shù)的取值,灰色區(qū)域是分布的置信區(qū)間(均值±標(biāo)準(zhǔn)差),點(diǎn)代表模式參數(shù)A的參考真值,虛線代表的是分布的均值；(a)代表由觀測(cè)區(qū)域a中的觀測(cè)反演的結(jié)果；(b)代表由觀測(cè)區(qū)域b中的觀測(cè)反演的結(jié)果Fig.10.The result retrieved by different observations.The horizontal coordinate represents the height,the vertical coordinate is the value of model parameter.The grey shaded area is the confidence interval of A(mean±standard deviation),the solid point is the reference of A,the dotted line is the mean of distribution.(a)display the result retrieved by observations in area a,(b)display the result retrieved by observations in area b.

5 總結(jié)與討論

本文利用PC-EnKF方法以及風(fēng)速的觀測(cè),對(duì)SGEM模式中湍流黏性系數(shù)的不確定性進(jìn)行了反演,求出了其后驗(yàn)概率密度函數(shù).采用兩類不同類型的湍流黏性系數(shù)進(jìn)行數(shù)值實(shí)驗(yàn),實(shí)驗(yàn)一中假定待反演的系數(shù)是常數(shù),實(shí)驗(yàn)二中假定是符合先驗(yàn)協(xié)方差函數(shù)的隨高度變化的系數(shù),均采用PC-EnKF進(jìn)行反演.實(shí)驗(yàn)結(jié)果表明,PC-EnKF能夠很有效地求解后驗(yàn)概率密度函數(shù),混沌多項(xiàng)式展開能夠有效進(jìn)行不確定性的表達(dá),將PC與EnKF結(jié)合起來能夠用于不確定性的反問題.而且我們?cè)趯?shí)驗(yàn)中發(fā)現(xiàn),根據(jù)系數(shù)的先驗(yàn)分布計(jì)算出風(fēng)速的先驗(yàn)分布,在此基礎(chǔ)上找到風(fēng)速不確定性大的區(qū)域,在此區(qū)域內(nèi)的觀測(cè)數(shù)據(jù)在反演系數(shù)時(shí)具有更好的效果,這對(duì)于觀測(cè)點(diǎn)的位置的選擇具有指導(dǎo)意義.

本文利用混沌多項(xiàng)式展開的方法代替EnKF中的蒙特卡羅法,能夠很有效地進(jìn)行不確定性問題的反問題.但是由于混沌多項(xiàng)式展開對(duì)于高維的不確定性問題的計(jì)算量急劇增加,所以PC-EnKF對(duì)于高維問題的求解還存在著一些問題,下一步的工作中將重點(diǎn)關(guān)注不確定性降維方面.