亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        HBFNEnKF混合同化方法設(shè)計(jì)及檢驗(yàn)

        2016-10-13 17:30:48朱浩楠閔錦忠杜寧珠
        大氣科學(xué) 2016年5期
        關(guān)鍵詞:效果方法

        朱浩楠 閔錦忠, 2 杜寧珠, 2

        ?

        HBFNEnKF混合同化方法設(shè)計(jì)及檢驗(yàn)

        朱浩楠1閔錦忠1, 2杜寧珠1, 2

        1南京信息工程大學(xué)氣象災(zāi)害預(yù)報(bào)預(yù)警與評(píng)估協(xié)同創(chuàng)新中心,南京210044;2南京信息工程大學(xué)氣象災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室,南京210044

        基于前后張馳逼近(Back and Forth Nudging,簡(jiǎn)稱BFN)和集合卡爾曼濾波(EnKF)方法,構(gòu)建了一種新的同化方法HBFNEnKF(Hybrid Back and Forth Nudging EnKF)混合同化方法,并將此同化系統(tǒng)分別與通道淺水模式(shallow water model)和全球淺水模式對(duì)接,檢驗(yàn)了HBFNEnKF同化方法的有效性。同時(shí),對(duì)比了集合均方根濾波(EnSRF)、HNEnKF (Hybrid Nudging EnKF)、HBFNEnKF三種方法在有誤差模式中的同化效果。試驗(yàn)結(jié)果表明:HBFNEnKF同化方法保留了HNEnKF方法的同化連續(xù)性,解決了EnKF同化不連續(xù)不平滑的問題,同時(shí)還有著更快的收斂速度;當(dāng)采用單變量分析試驗(yàn)時(shí),HBFNEnKF方法的優(yōu)勢(shì)最為明顯,表明HBFNEnKF能夠較好地保持不同模式變量間的平衡。此外,增量場(chǎng)尺度分析結(jié)果表明:相比EnSRF,HBFNEnKF在大尺度范圍有更好的同化效果,同時(shí)能夠避免在中小尺度范圍內(nèi)出現(xiàn)大量的虛假增量。

        資料同化 Hybrid Back and Forth Nudging EnKF (HBFNEnKF) Ensemble Square-Root Fiter (EnSRF) 淺水模式

        1 引言

        數(shù)值天氣預(yù)報(bào)是一個(gè)初始場(chǎng)和邊界條件的問題。準(zhǔn)確的初始場(chǎng)在提高數(shù)值天氣預(yù)報(bào)中起著重要作用,資料同化是為數(shù)值模式提供接近實(shí)際大氣狀態(tài)初始場(chǎng)的重要手段,其主要目的是通過適宜的算法將新的觀測(cè)資料融入到模式場(chǎng)中,以得到更好的模式初值(Lewis et al., 2006)。經(jīng)過數(shù)十年的發(fā)展,當(dāng)前主流的同化方法分為變分同化法和集合同化法。前者主要包含三維變分(3DVAR)和四維變分(4DVAR)(Lewis et al., 2006);后者主要包含以集合卡爾曼濾波(EnKF)(Evensen, 1994, 2003)為核心的,經(jīng)過不同拓展延伸的集合同化方法,如集合變換卡爾曼濾波(ETKF)(Bishop et al., 2001; Livings, 2005)、集合調(diào)整卡爾曼濾波(EAKF)(Anderson, 2001)、集合均方根濾波(EnSRF)(Whitaker and Hamill, 2002; 邵長(zhǎng)亮和閔錦忠,2015)以及很多其他改進(jìn)算法如迭代集合平方根濾波(iEnSRF)(閔錦忠等,2012;王世璋等,2013)等。

        盡管同化方法眾多,但每一種同化思想都有各自的優(yōu)缺點(diǎn)。EnKF方法由于采用了集合預(yù)報(bào)的思想,利用集合成員來計(jì)算背景誤差協(xié)方差矩陣(矩陣),實(shí)現(xiàn)了背景誤差協(xié)方差的流依賴,同時(shí)避免了4DVAR同化中對(duì)切線性模式和伴隨模式的構(gòu)建。這使得EnKF相比3DVAR、4DVAR更加適合應(yīng)用到發(fā)展較快、變化迅速的中小尺度過程中(Yang et al., 2009)。然而,為了構(gòu)造具有代表性的矩陣需要大量的集合成員,增大了模式的計(jì)算量;而較小的集合成員容易產(chǎn)生濾波發(fā)散。同時(shí),EnKF方法的局地化操作以及在同化過程中造成的不連續(xù)性所帶來的高頻振蕩可能會(huì)破壞模式變量的平衡關(guān)系(Bloom et al., 1996; Ourmières et al., 2006; Kepert, 2009)并造成數(shù)據(jù)損失,從而引入誤差。這些缺點(diǎn)在對(duì)初值敏感的中小尺度過程中都是極為不利的。

        張馳逼近(Nudging)(Hoke and Anthes, 1976)是一種傳統(tǒng)的連續(xù)同化方法,通過對(duì)模式積分方程添加Nudging項(xiàng)使每一步積分都能向觀測(cè)場(chǎng)逼近,并且積分時(shí)間越長(zhǎng)同化效果越好(Lakshmivarahan and Lewis, 2013)。Nudging方法作用于每一積分步,其同化過程不會(huì)造成嚴(yán)重的不連續(xù)問題;由于利用了模式積分進(jìn)行動(dòng)力約束,不會(huì)破壞模式各變量間的物理平衡;但存在同化收斂較慢、資料利用不夠充分的問題。為了克服這一缺點(diǎn),Auroux and Blum(2005, 2008)提出了Back and Forth Nudging(BFN)方法:通過向前、向后積分迭代增加積分時(shí)間,以加快Nudging同化收斂速度,實(shí)現(xiàn)對(duì)觀測(cè)資料的充分利用。通過與中尺度氣象業(yè)務(wù)預(yù)報(bào)模式Meso-NH model(Boilley and Mahfouf, 2012)和海洋業(yè)務(wù)預(yù)報(bào)模式NEMO(Ruggiero et al., 2015)對(duì)接,試驗(yàn)證明BFN方法有與4DVAR相當(dāng)?shù)耐Ч?/p>

        為了解決EnKF的同化不連續(xù)問題,Lei et al.(2012a, 2012b)提出了Hybrid Nudging EnKF(HNEnKF)方法,其基本思想是使用EnKF的卡爾曼增益矩陣作為Nudging算子引入模式預(yù)報(bào)方程,最后將Nudging同化的結(jié)果與利用EnKF更新后的集合擾動(dòng)相加,得到新的集合成員。Lei et al.(2012a, 2012b)在文章中用Lorenz 63模式和淺水模式進(jìn)行了試驗(yàn),取得了較好的效果:在Lorenz63模式中,HNEnKF方法的連續(xù)性好于傳統(tǒng)EnKF方法和分析增量法(IAU)(Bloom et al., 1996; Ourmières et al., 2006),并與集合卡爾曼平滑(EnKS)(Evensen and van Leeuwen, 2000)相當(dāng);在淺水模式下,HNEnKF的同化效果優(yōu)于EnKF和EnKS。然而當(dāng)同化窗口增大時(shí),HNEnKF方法 的同化效果明顯降低;模式均方根誤差(RMSE)收斂速度依然較慢;如本文前述,Nudging算法 的特性在于積分時(shí)間越長(zhǎng),同化效果越好(Lakshmivarahan and Lewis, 2013),簡(jiǎn)單地將Nudging算子用卡爾曼增益矩陣進(jìn)行替換的混合方法,并不能充分發(fā)揮其同化效果。

        為解決上述問題,本文根據(jù)HNEnKF方法的混合思路,將BFN算法與EnKF進(jìn)行結(jié)合,設(shè)計(jì)了新的混合算法:Hybrid Back and Forth Nudging EnKF(HBFNEnKF)同化方法,并用淺水模式(shallow water model)進(jìn)行了幾組同化試驗(yàn),詳細(xì)討論了結(jié)果。新算法采用了BFN前后積分迭代的思路,使同化過程能很快收斂,同時(shí)能對(duì)非觀測(cè)模式變量進(jìn)行較好調(diào)整;此外,HBFNEnKF保留了HNEnKF良好的同化連續(xù)性;在模式方程約束下,其各變量間的平衡關(guān)系得到維持。對(duì)HNEnKF和HBFNEnKF方法的介紹將在文章第二部分給出,第三部分的內(nèi)容是試驗(yàn)方案設(shè)計(jì),試驗(yàn)結(jié)果與分析在第四部分給出,第五部分為全文總結(jié)。

        2 基本同化算法介紹

        2.1 EnSRF同化方法

        EnSRF是EnKF同化方法的一個(gè)重要分支,本文設(shè)計(jì)的HBFNEnKF方法也采用了EnSRF的思路計(jì)算卡爾曼增益矩陣、更新集合擾動(dòng)。傳統(tǒng)的EnKF為了避免對(duì)背景誤差協(xié)方差的低估,采用了在觀測(cè)資料中添加隨機(jī)擾動(dòng)的方法,這樣一來會(huì)在同化中引入擾動(dòng)誤差(Whitaker and Hamill, 2002)。EnSRF方案則通過一系列的變換,避免了這一誤差。

        EnSRF的同化過程分為兩步:

        其中,(1)式為預(yù)報(bào)步,代表模式預(yù)報(bào)方程,代表模式預(yù)報(bào)變量,上標(biāo)a代表同化后的分析值,上標(biāo)b代表預(yù)報(bào)值或背景場(chǎng),下標(biāo)代表第個(gè)集合成員。(2)、(3)式為分析步,(2)式用于更新集合平均,(3)式用于更新集合擾動(dòng),為觀測(cè)算子,代表觀測(cè)場(chǎng),是EnSRF為了避免EnKF中對(duì)觀測(cè)場(chǎng)加擾而引入的參數(shù),為卡爾曼增益矩陣。

        2.2 Nudging同化方法和BFN同化方法介紹

        Nudging同化也被稱為牛頓張馳逼近法,是一種連續(xù)同化方法,通過在預(yù)報(bào)方程加入Nudging項(xiàng)使得模式變量在積分過程中逐步逼近真值。傳統(tǒng)的Nudging同化公式為

        其中,

        圖1 權(quán)重系數(shù)wt隨時(shí)間變化的示意圖

        理論上,采用Nudging同化會(huì)使模式場(chǎng)隨著積分不斷收斂于真值(Lakshmivarahan and Lewis, 2013),但要求Nudging算子既不能太大,導(dǎo)致積分不穩(wěn)定,也不能太小,使得收斂速度太慢,降低同化效果。在通常的業(yè)務(wù)模式如WRF模式中,的取法為Cressman反距離權(quán)重法(Stauffer and Seaman, 1990; Liu et al., 2005, 2006),Stauffer and Seaman(1990),Stauffer and Bao(1993)以及Zou et al.(1992)提出了最優(yōu)逼近(Optimal Nudging),利用伴隨矩陣求解對(duì)應(yīng)的代價(jià)函數(shù),得出某一同化區(qū)間的最優(yōu)Nudging算子。然而這樣的求解方式需要伴隨模式的構(gòu)建,加大了操作難度。

        從增加模式積分時(shí)間的角度考慮,Auroux and Blum(2005, 2008)提出了BFN方法,其基本思路為

        (11)式為向前逼近(Forward Nudging),(12)式為向后逼近(Backward Nudging)。其中,、分別表示向前積分和向后積分變量,下標(biāo)表示迭代的次數(shù),、分別為向前Nudging算子和向后Nudging算子。Auroux and Blum(2005, 2008)將Nudging算子、處理為同一形式:

        BFN同化方法的步驟為:利用Forward Nudging將模式向前積分到同化窗口終點(diǎn),將模擬結(jié)果賦給Backward Nudging作初值;利用Backward Nudging將模式向后積分到窗口起始時(shí)刻,把模擬結(jié)果賦給Forward Nudging,由此不斷迭代直至收斂。由于模式的非線性特性,向后積分會(huì)產(chǎn)生較強(qiáng)的不穩(wěn)定,但在Backward Nudging里,Nudging項(xiàng)的引入使向后積分得以較為穩(wěn)定地進(jìn)行,由此構(gòu)建出不需要伴隨模式的四維同化方案。

        Auroux and Blum(2005, 2008)證明了BFN同化方法的收斂性和較好的同化效果,但實(shí)際上該方法也存在很多問題:(1)Auroux and Blum(2005, 2008)給出的Nudging算子僅是一個(gè)近似,并未給出明確的標(biāo)準(zhǔn),對(duì)的取值較為隨意;(2)觀測(cè)誤差協(xié)方差矩陣在實(shí)際中有時(shí)無法確定,而計(jì)算切線性觀測(cè)算子需要復(fù)雜的操作;(3)為了達(dá)到好的同化效果,BFN同化需要較多次迭代才能收斂(Auroux, 2009),使同化花費(fèi)了很多時(shí)間。

        2.3 混合方法介紹

        為了解決EnKF的同化不連續(xù)問題,Lei et al.(2012a, 2012b)提出了HNEnKF同化方案,該方法中Nudging算子的替換公式為

        HNEnKF的同化步驟為:首先將集合成員積分到同化時(shí)次,利用(4)式計(jì)算得到,然后按照(14)式的混合方案計(jì)算出Nudging對(duì)應(yīng)的同化算子,將模式積分到同化時(shí)刻;用Nudging積分后的結(jié)果替換掉集合平均場(chǎng),與EnKF更新后的集合擾動(dòng)相加,得到新的集合成員;將集合成員繼續(xù)向前積分,循環(huán)同化。HNEnKF與IAU有很多相似之處,但其最大不同在于:IAU在積分步上疊加的增量是靜態(tài)的,HNEnKF同化則使用隨時(shí)間變化的分析增量進(jìn)行疊加。圖2給出了HNEnKF方法的詳細(xì)同化流程。

        圖2 HNEnKF方案的同化流程示意圖

        如前所述,HNEnKF方法還存在一些問題。為了改進(jìn)該方法,本文結(jié)合Auroux and Blum(2005, 2008)提出的BFN算法與EnKF,設(shè)計(jì)了HBFNEnKF混合方法。其原理為:將(11)、(12)式中的、利用下式進(jìn)行替換:

        圖3 HBFNEnKF方案的同化流程示意圖

        3 試驗(yàn)設(shè)計(jì)

        為了檢驗(yàn)HBFNEnKF的同化效果,本文首先在通道淺水模式下進(jìn)行了同化試驗(yàn)。模式方程組為

        其中,、、的取值范圍滿足、、,和表示模式緯向和經(jīng)向距離,為模式積分時(shí)間,、表示水平風(fēng)場(chǎng),和s分別代表模式高度場(chǎng)和模式地形,代表粘性耗散系數(shù),、分別代表每個(gè)模式變量對(duì)應(yīng)的Nudging算子。模式考慮平面假設(shè),科氏參數(shù)取為常數(shù)10?4s?1,重力加速度取為9.8 m s?1,格距?=?=30 km,=1500 km,=1500 km,=105m2s?1。模式格點(diǎn)采用荒川c網(wǎng)格,積分步長(zhǎng)?= 360 s,積分方法采用龍格庫塔4階積分方法。模擬真實(shí)場(chǎng)時(shí)地形設(shè)置為

        而為了模擬模式誤差,在同化試驗(yàn)中將地形設(shè)置為

        (19)式中,、、的取值范圍滿足、、,(20)式中,、、的取值范圍滿足、、。、、代表向前積分的模式變量,、、代表向后積分的模式變量,下標(biāo)代表迭代次數(shù),、、分別代表不同模式變量對(duì)應(yīng)的向前積分Nudging算子,、、分別代表不同模式變量對(duì)應(yīng)的向后積分Nudging算子。

        模式東西邊界采用周期邊界條件,上下邊界設(shè)為0,初始高度場(chǎng)利用下式得到:

        其中,、滿足的取值范圍為、,=3000 m,,和分別代表緯向距離和經(jīng)向距離。計(jì)算得到高度場(chǎng)后,再利用地轉(zhuǎn)關(guān)系計(jì)算出初始風(fēng)場(chǎng)。

        同化試驗(yàn)使用了50個(gè)集合,其構(gòu)造方式為:在初始高度場(chǎng)上疊加滿足分布的擾動(dòng),然后利用地轉(zhuǎn)關(guān)系計(jì)算出對(duì)應(yīng)風(fēng)場(chǎng)集合成員。真實(shí)場(chǎng)和集合成員都經(jīng)過60 h的啟動(dòng)時(shí)間(spin up),同時(shí)也是對(duì)集合成員的發(fā)展。同化時(shí)間間隔分為24 h一次和48 h一次兩種。生成含誤差觀測(cè)資料的方法為:將高度場(chǎng)和風(fēng)場(chǎng)、三個(gè)變量的真值插值到模式空間中120個(gè)不規(guī)則分布點(diǎn)上,再分別疊加滿足分布(高度場(chǎng))的擾動(dòng)和滿足分布(風(fēng)場(chǎng))的擾動(dòng)得到。

        試驗(yàn)中所用同化方案為EnSRF、HNEnKF以及HBFNEnKF,并另外進(jìn)行一次模式自由積分作為控制試驗(yàn)(CTRL)。所有同化試驗(yàn)均采用GC(Gaspari and Cohn, 1999)局地化方法進(jìn)行局地化處理,局地化半徑取為600 km;同時(shí)采用松弛膨脹法進(jìn)行協(xié)方差膨脹(relax inflation)(Zhang et al., 2004),膨脹系數(shù)設(shè)置為0.1。利用本文前述的時(shí)間權(quán)重算子,將HBFNEnKF和HNEnKF中的Nudging項(xiàng)分配到以觀測(cè)資料加入時(shí)次為中心的、左右各1 h的時(shí)間窗口內(nèi)。

        為了檢驗(yàn)各方案同化過程中變量間的平衡性,上述試驗(yàn)均關(guān)閉了EnSRF同化中的多變量分析。所謂多變量分析,指集合卡爾曼濾波通過背景誤差協(xié)方差,在同化時(shí)實(shí)現(xiàn)各模式變量相互影響;關(guān)閉多變量分析,即進(jìn)行單一變量要素同化試驗(yàn)。為了進(jìn)一步檢驗(yàn)新方案的效果,本文還在全球淺水模式中進(jìn)行了考慮多變量分析后的同化試驗(yàn)。全球淺水模式緯向格點(diǎn)180個(gè),經(jīng)向格點(diǎn)90個(gè),積分步長(zhǎng)60 s,模式初始場(chǎng)設(shè)置參見Nair et al.(2005)。觀測(cè)資料取為405個(gè)不規(guī)則分布的點(diǎn)。各試驗(yàn)組都取20個(gè)集合成員,局地化半徑為50 km,協(xié)方差膨脹系數(shù)為0.3。在全球淺水模式試驗(yàn)中,首先讓模式積分10 d作為spin up并發(fā)展集合成員,再進(jìn)行3 d、每6 h一次的同化循環(huán),之后再進(jìn)行4 d的自由積分。

        所有試驗(yàn)方案詳細(xì)配置參見表1和表2。

        表1 通道淺水模式中各試驗(yàn)組名稱及設(shè)置

        Table 1 Names and designs of experiments in the channel shallow water model

        表2 全球淺水模式中各試驗(yàn)組名稱及設(shè)置

        Table 2 Names and designs of experiments in the global shallow water model

        4 試驗(yàn)結(jié)果與分析

        4.1 通道淺水模式試驗(yàn)結(jié)果和分析

        圖4是同化間隔為24 h的高度場(chǎng)、風(fēng)場(chǎng)均方根誤差(RMSE)隨時(shí)間的變化趨勢(shì)。圖4a、b、c分別為僅同化高度場(chǎng)、僅同化風(fēng)場(chǎng)、同時(shí)同化高度場(chǎng)和風(fēng)場(chǎng)時(shí)的高度場(chǎng)RMSE變化;圖4d、e、f為對(duì)應(yīng)的風(fēng)場(chǎng)RMSE變化。由于與Nudging混合,HNEnKF和HBFNEnKF方法單次同化增量很小,這里給出的是每一積分步的RMSE??梢悦黠@看出,不論是同化單一變量觀測(cè)還是同時(shí)同化高度場(chǎng)和風(fēng)場(chǎng)觀測(cè),HBFNEnKF的同化效果都是最好的。

        圖4 通道淺水模式下各同化方案高度場(chǎng)的RMSE隨時(shí)間的變化:(a)只同化高度場(chǎng);(b)只同化風(fēng)場(chǎng);(c)同時(shí)同化高度場(chǎng)和風(fēng)場(chǎng)。(d–f)同(a–c),但為風(fēng)場(chǎng)的RMSE變化

        當(dāng)僅同化高度場(chǎng)時(shí)(圖4a),各方案的同化效果都較差,但EnSRF在每個(gè)同化時(shí)刻后都有很強(qiáng)的波動(dòng),其RMSE有先增大再逐漸減小的變化過程;同樣的現(xiàn)象也出現(xiàn)在圖4b中。這是由于只采用了單變量分析,EnSRF同化不能很好地使模式各變量之間達(dá)到平衡關(guān)系,需要一段時(shí)間的spin up。該問題并未出現(xiàn)在HNEnKF和HBFNEnKF試驗(yàn)中,不論是只同化高度場(chǎng)資料還是只同化風(fēng)場(chǎng)資料,兩種混合方法都能保持較好的同化連續(xù)性,并維持模式各變量間的平衡關(guān)系。由圖4a和圖4d可見,由于只同化了高度場(chǎng)資料,HNEnKF雖然比EnSRF的同化效果好,但其優(yōu)勢(shì)并不明顯,且需要在積分一段時(shí)間后,才能加大與EnSRF同化間的差距;而HBFNEnKF不但能迅速訂正高度場(chǎng),同時(shí)能很好地改善風(fēng)場(chǎng)。這說明經(jīng)過前后積分迭代,HBFNEnKF比HNEnKF能更好地進(jìn)行觀測(cè)資料和模式變量間的調(diào)整。

        圖5為模式循環(huán)同化期間,各同化方法在不同同化時(shí)間間隔情況下時(shí)間平均的RMSE和不連續(xù)性參數(shù)d(discontinue parameter)的對(duì)比。這里d為L(zhǎng)ei et al.(2012a, 2012b)引入的參數(shù),目的是為了檢驗(yàn)?zāi)骋煌椒ǖ耐B續(xù)性。其計(jì)算公式為

        其中,為試驗(yàn)中同化時(shí)刻的總次數(shù),E-1和E分別表示某同化時(shí)刻同化前和同化后的RMSE。d值越大則說明連續(xù)性越差。

        圖5 模式循環(huán)同化期間平均的(a)高度場(chǎng)、(b)風(fēng)場(chǎng)的RMSE和(c)高度場(chǎng)、(d)風(fēng)場(chǎng)的Pd分布示意圖。(c)、(d)中,HNEnKF和HBFNEnKF試驗(yàn)的值在右側(cè)坐標(biāo)中顯示

        由圖5a、b可見,HBFNEnKF試驗(yàn)的時(shí)間平均RMSE最小,HNEnKF試驗(yàn)次之,EnSRF試驗(yàn)最大。隨著同化時(shí)間間隔的增大,各同化方法的RMSE都有增加,其中變化最明顯的是HNEnKF方法,而同樣條件下HBFNEnKF方案的RMSE增幅并不明顯。根據(jù)同化間隔為48 h時(shí)的RMSE變化曲線(圖略)可知,相比間隔時(shí)間為24 h的情況,各同化方法的均方根誤差收斂時(shí)間都有增加,這是時(shí)間平均RMSE增加的主要原因,但其中以HBFNEnKF方法的收斂時(shí)間變化最小。此外,HNEnKF與EnSRF之間的RMSE差距明顯減小。Lei et al.(2012a, 2012b)指出這是因?yàn)闀r(shí)間窗口增大,Nudging同化受到更多模式誤差帶來的影響。而采用前后積分迭代進(jìn)行多次調(diào)整,HBFNEnKF在一定程度上改善了HNEnKF這一缺點(diǎn)。由圖5c、d可見,HBFNEnKF同化也保留了HNEnKF良好的同化連續(xù)性;與他們相比,EnSRF的同化連續(xù)性最差。

        圖6為同化時(shí)間間隔24 h,只同化風(fēng)場(chǎng)資料時(shí),各試驗(yàn)循環(huán)同化結(jié)束后的環(huán)流形勢(shì)??梢钥闯?,循環(huán)同化結(jié)束后模式場(chǎng)中主要的系統(tǒng)是位于=15 km處的兩個(gè)高壓中心以及位于=30 km處的兩個(gè)小低壓。由圖6c、d、e可見,EnSRF試驗(yàn)未能將=30 km處的小低壓模擬出來。同時(shí)兩個(gè)高壓中心的強(qiáng)度、位置以及等高線的流形都與真實(shí)場(chǎng)存在一定差異。雖然沒有模擬出低壓中心,但是在真實(shí)場(chǎng)中小低壓的位置(=30 km,=20 km),EnSRF模擬出了氣旋性環(huán)流。這說明EnSRF同化風(fēng)場(chǎng)資料時(shí),不能很好地調(diào)整模式高度場(chǎng)。由圖6d可見,相比于EnSRF,HNEnKF有一定改進(jìn),但依然不足;同化效果最好的是HBFNEnKF試驗(yàn),從各環(huán)流系統(tǒng)的強(qiáng)度、位置分布看,HBFNEnKF都有與真實(shí)場(chǎng)最接近的結(jié)果。

        圖6 只同化風(fēng)場(chǎng)、同化周期24 h時(shí),各試驗(yàn)循環(huán)同化結(jié)束后高度場(chǎng)(單位:m)和風(fēng)場(chǎng)(單位:m s?1):(a)真實(shí)場(chǎng);(b)控制試驗(yàn);(c)EnSRF;(c)HNEnKF;(e)HBFNEnKF

        4.2 全球淺水模式試驗(yàn)結(jié)果和分析

        圖7給出的是考慮了多變量分析,各同化方案在全球淺水模式同化試驗(yàn)中RMSE的變化趨勢(shì)。由于全球淺水模式試驗(yàn)積分步較多,所以制作插圖時(shí)只考慮了同化時(shí)刻。可以明顯發(fā)現(xiàn),在考慮多變量分析后,EnSRF方案的同化效果得到了很大改進(jìn)。然而其同化不連續(xù)、不平衡的現(xiàn)象依然存在。在3 d的循環(huán)同化(圖中黑色豎線)后,EnSRF試驗(yàn)的RMSE有較快增加;而從整體來看,HNEnKF和HBFNEnKF在自由積分階段的RMSE都比EnSRF低。從整體來看,HBFNEnKF依然保留了收斂快、同化效果好、RMSE低的特點(diǎn)。不論是同化單一觀測(cè)變量還是同時(shí)同化風(fēng)場(chǎng)和高度場(chǎng)變量,HBFNEnKF都能保持最小的均方根誤差。尤其是在只同化風(fēng)場(chǎng)觀測(cè)和只同化高度場(chǎng)觀測(cè)時(shí),HBFNEnKF對(duì)非觀測(cè)模式變量能有較好訂正(圖7b、d)。

        圖7 全球淺水模式下各同化方案高度場(chǎng)的RMSE隨時(shí)間的變化:(a)只同化高度場(chǎng);(b)只同化風(fēng)場(chǎng)、(c)同時(shí)同化風(fēng)場(chǎng)和高度場(chǎng)觀測(cè)。(d–f)同(a–c),但為風(fēng)場(chǎng)的RMSE變化

        均方根誤差只是對(duì)同化效果整體情況的反映。圖8、9給出了在全球淺水模式中各試驗(yàn)循環(huán)同化結(jié)束時(shí),且自由積分4 d后的環(huán)流形勢(shì)(180°為圖的中心經(jīng)度)。圖8是只同化風(fēng)場(chǎng)觀測(cè)的結(jié)果;圖9為只同化高度場(chǎng)觀測(cè)的結(jié)果。

        圖8 全球淺水模式下只同化風(fēng)場(chǎng)資料時(shí),各同化方案循環(huán)同化結(jié)束后,并自由積分4天后的環(huán)流形勢(shì):(a)真實(shí)場(chǎng);(b)EnSRF;(c)HNEnKF;(d)HBFNEnKF。彩色陰影代表位勢(shì)高度(單位:m),黑色線代表等高線(單位:m),箭頭代表風(fēng)場(chǎng)(單位:m s?1)

        圖9 同圖8,但為只同化高度場(chǎng)資料的結(jié)果

        由圖8可見,當(dāng)只同化風(fēng)場(chǎng)資料時(shí),HBFNEnKF相對(duì)于另外兩個(gè)試驗(yàn),在各環(huán)流系統(tǒng)中心的強(qiáng)度和位置上都與真實(shí)場(chǎng)最為接近,較好地模擬出了位于30°~60°N的一高一低兩個(gè)環(huán)流中心,以及位于30°~60°S、120°W處的高壓中心。HNEnKF和EnSRF試驗(yàn)在30°N、150°W處模擬出了一個(gè)虛假的高壓中心,而HBFNEnKF試驗(yàn)的模擬結(jié)果則沒有出現(xiàn)這樣的虛假高壓中心。就風(fēng)場(chǎng)的分布和強(qiáng)度來說,HBFNEnKF的同化效果也優(yōu)于另外兩種方案。

        由圖9可見,當(dāng)只同化高度場(chǎng)資料時(shí),EnSRF和HNEnKF試驗(yàn)?zāi)M結(jié)果中的虛假環(huán)流系統(tǒng)發(fā)展得更為明顯。相比只同化風(fēng)場(chǎng)的試驗(yàn),由于觀測(cè)資料和觀測(cè)變量的減少,HBFNEnKF同化效果有一定程度的降低。盡管如此,試驗(yàn)HBFNEnKF模擬的主要環(huán)流系統(tǒng)及其強(qiáng)度與位置依然優(yōu)于另外兩組試驗(yàn)。

        4.3 增量場(chǎng)分析

        本節(jié)將通過增量場(chǎng)的分析簡(jiǎn)單說明,相比EnSRF,HBFNEnKF能夠取得更好的同化效果,同時(shí)還能避免一些虛假環(huán)流系統(tǒng)的出現(xiàn)。由于HBFNEnKF能夠通過反向積分回到初始時(shí)刻,求得對(duì)應(yīng)時(shí)刻同化前后的增量,而HNEnKF僅為向前積分且單次同化增量過小,無法進(jìn)行比較,所以這里只用HBFNEnKF與EnSRF進(jìn)行對(duì)比。

        首先利用離散余弦變化(DCT)(Denis et al., 2002)將EnSRF、HBFNEnKF以及真實(shí)場(chǎng)在模式初始時(shí)刻的高度場(chǎng)轉(zhuǎn)換到波數(shù)空間,然后利用帶通濾波對(duì)其進(jìn)行尺度分離,最后求得各尺度上的增量場(chǎng)。圖10為第一個(gè)同化時(shí)刻、同時(shí)同化高度場(chǎng)和風(fēng)場(chǎng)時(shí),經(jīng)過尺度分離后,EnSRF和HBFNEnKF高度場(chǎng)的增量場(chǎng)以及真實(shí)場(chǎng)與相應(yīng)時(shí)刻控制試驗(yàn)的差異場(chǎng)(因?yàn)槭堑谝粋€(gè)同化時(shí)刻,所以EnSRF和HBFNEnKF同化前的背景場(chǎng)與相應(yīng)時(shí)刻的控制試驗(yàn)是相同的),在小尺度、中尺度以及大尺度范圍內(nèi)的差異分布。可以發(fā)現(xiàn):在大尺度范圍(圖10c、f、i),不論是增量的強(qiáng)度還是分布位置,HBFNEnKF都比EnSRF更接近真實(shí)場(chǎng);在中尺度范圍(圖10b、e、h),雖然HBFNEnKF的增量場(chǎng)與真實(shí)場(chǎng)和控制試驗(yàn)的差異場(chǎng)沒有很好的匹配,但是與EnSRF的增量場(chǎng)相比,避免了許多虛假增量的出現(xiàn);在小尺度范圍(圖10a、d、g),HBFNEnKF的增量場(chǎng)與真實(shí)場(chǎng)和控制試驗(yàn)的差異場(chǎng)分布較為一致,EnSRF有明顯的虛假增量,這正是由于EnSRF同化不平衡、不連續(xù)所致。

        圖10 全球淺水模式下第一個(gè)同化時(shí)刻同時(shí)同化高度場(chǎng)和風(fēng)場(chǎng)時(shí),(a)EnSRF、(d)HBFNEnKF、(g)真實(shí)場(chǎng)與控制試驗(yàn)在小尺度范圍上高度場(chǎng)的差異。(b、e、h)同(a、d、g),但為在中尺度范圍上的差異場(chǎng);(c、f、i)同(a、d、g),但為在大尺度范圍上的差異場(chǎng)

        圖11為第一個(gè)同化時(shí)刻、只同化高度場(chǎng)時(shí),進(jìn)行了尺度分離后,EnSRF、HBFNEnKF以及真實(shí)場(chǎng)的緯向風(fēng)分量與控制試驗(yàn)的差異場(chǎng)。只同化高度場(chǎng)是為了使差異更加明顯,同時(shí)體現(xiàn)同化方案對(duì)非觀測(cè)模式變量的訂正效果??梢园l(fā)現(xiàn):在大尺度范圍,HBFNEnKF有著與真實(shí)場(chǎng)更為一致的增量分布;而在中尺度尤其是小尺度范圍,HBFNEnKF避免了在EnSRF同化中出現(xiàn)的虛假增量。這說明,在對(duì)非觀測(cè)模式變量的訂正上,HBFNEnKF有比EnSRF更優(yōu)秀的能力,同時(shí)能解決EnSRF同化中的虛假相關(guān)問題。

        圖11 同圖10,但為只同化高度場(chǎng)時(shí)緯向風(fēng)分量的結(jié)果

        綜上所述,HBFNEnKF同化能在模式大尺度范圍進(jìn)行更好的同化,同時(shí)在中小尺度范圍能夠避免傳統(tǒng)EnSRF同化中虛假相關(guān)的出現(xiàn)。這也正是HBFNEnKF有比EnSRF更好的同化效果、并能夠避免模式結(jié)果中虛假環(huán)流系統(tǒng)出現(xiàn)的原因。

        5 結(jié)論與討論

        本文基于Lei et al.(2012a, 2012b)設(shè)計(jì)的HNEnKF同化方法和Auroux and Blum(2005, 2008)設(shè)計(jì)的BFN方法,通過結(jié)合BFN和EnKF同化算法,設(shè)計(jì)了HBFNEnKF混合同化方案,并針對(duì)EnSRF、HNEnKF和HBFNEnKF三種算法,在考慮了誤差的通道淺水模式、全球淺水模式中進(jìn)行同化試驗(yàn)。試驗(yàn)對(duì)比了在單變量分析、多變量分析情況下,這三種同化方法的RMSE變化、不連續(xù)性參數(shù)分布以及經(jīng)過循環(huán)同化并自由積分后的模式環(huán)流形勢(shì),得出以下結(jié)論:

        (1)HBFNEnKF混合同化方案保留了HNEnKF同化的連續(xù)性和平滑性,解決了EnSRF同化不連續(xù)、不平滑的缺點(diǎn)。

        (2)通過前后迭代積分,HBFNEnKF充分發(fā)揮了Nudging同化的特點(diǎn):積分時(shí)間越長(zhǎng)同化效果越好,使得同化能很快收斂,相比EnSRF和HNEnKF有最快的收斂速度,同時(shí)取得了更好的同化效果。而通過前后積分迭代,HBFNEnKF能較好地逼近真實(shí)場(chǎng),并在一定程度上改善了HNEnKF同化效果隨著同化時(shí)間間隔增加而降低的問題。

        (3)由于HBFNEnKF作用于模式積分,所以HBFNEnKF同化能利用模式方程進(jìn)行約束,使得變量間的平衡得以維持。在單變量分析試驗(yàn)中,HBFNEnKF和HNEnKF都能較好地利用觀測(cè)資料對(duì)各模式變量進(jìn)行訂正。而相比之下,由于沒有考慮多變量分析,EnSRF方法在同化單一變量時(shí)帶來了較大的模式擾動(dòng),增加了模式的spin up時(shí)間,從而影響了同化效果。

        (4)全球淺水模式的試驗(yàn)表明,在經(jīng)過4 d自由積分后,HBFNEnKF不但更加接近真實(shí)場(chǎng),同時(shí)也避免了一些虛假環(huán)流系統(tǒng)的出現(xiàn)。通過對(duì)高度場(chǎng)增量的尺度分析發(fā)現(xiàn),這是由于HBFNEnKF在大尺度范圍有比EnSRF更好的同化效果:其增量場(chǎng)在強(qiáng)度和空間上都更加接近同化前與真實(shí)場(chǎng)的差異場(chǎng)。在中小尺度范圍,HBFNEnKF還避免了一些虛假增量的出現(xiàn),這正是由于HBFNEnKF較好的同化連續(xù)性和平滑性所致。僅同化高度場(chǎng)時(shí),緯向風(fēng)增量的尺度分析也有一致結(jié)論,這是HBFNEnKF能夠較好地維持模式變量間相互平衡、對(duì)非觀測(cè)模式變量能進(jìn)行合理調(diào)整的證明。

        盡管HBFNEnKF在本次試驗(yàn)中有最佳的效果,但是淺水模式只是一個(gè)二維理想模式,并不能完全代表真實(shí)大氣的運(yùn)動(dòng)與結(jié)構(gòu)。HBFNEnKF同化方法在多層模式以及實(shí)際預(yù)報(bào)模式中的試驗(yàn)將在之后進(jìn)行,對(duì)真實(shí)觀測(cè)資料的同化效果也有待檢驗(yàn)。另外,HBFNEnKF的許多特性還有待進(jìn)一步探討,如:HBFNEnKF是否也保留了BFN能將大氣低層資料同化到模式高層的特性?迭代過程中的選取對(duì)同化效果有何影響?此外,Lorenz 96模式的試驗(yàn)結(jié)果(略)表明,通過引入不隨時(shí)間變化的部分,能夠進(jìn)一步增加HBFNEnKF的同化效果,這一改進(jìn)是否能在更為復(fù)雜的模式中實(shí)現(xiàn)?同時(shí),HBFNEnKF還存在由于迭代帶來的更多計(jì)算量,以及向后積分不穩(wěn)定等問題,雖然在本文試驗(yàn)中并未出現(xiàn)上述情況,但依然應(yīng)該進(jìn)行相應(yīng)的研究和改進(jìn)。

        參考文獻(xiàn)(:References)

        Anderson J L. 2001. An ensemble adjustment Kalman filter for data assimilation [J]. Mon. Wea. Rev., 129 (12): 2884–2903, doi:10.1175/ 1520-0493(2001)129<2884:AEAKFF>2.0.CO;2.

        Auroux D. 2009. The back and forth nudging algorithm applied to a shallow water model, comparison and hybridization with the 4D-VAR [J]. Int. J. Numer. Methods Fluids, 61 (8): 911–929, doi:10.1002/fld.1980.

        Auroux D, Blum J. 2005. Back and forth nudging algorithm for data assimilation problems [J]. Comptes Rendus Mathematique, 340 (12): 873–878, doi:10.1016/j.crma.2005.05.006.

        Auroux D, Blum J. 2008. A nudging-based data assimilation method: The Back and Forth Nudging (BFN) algorithm [J]. Nonlin. Processes Geophys., 15 (2): 305–319, doi:10.5194/npg-15-305-2008.

        Bishop C H, Etherton B J, Majumdar S J. 2001. Adaptive sampling with the ensemble transform Kalman filter. Part I: Theoretical aspects [J]. Mon. Wea. Rev., 129 (3): 420–436, doi:10.1175/1520-0493(2001)129<0420: ASWTET>2.0.CO;2.

        Bloom S C, Takacs L L, da Silva A M, et al. 1996. Data assimilation using incremental analysis updates [J]. Mon. Wea. Rev., 124 (6): 1256–1271, doi:10.1175/1520-0493(1996)124<1256:DAUIAU>2.0.CO;2.

        Boilley A, Mahfouf J F. 2012. Assimilation of low-level wind in a high- resolution mesoscale model using the back and forth nudging algorithm [J]. Tellus A, 64: 18697, doi:10.3402/tellusa.v64i0.18697.

        Denis B, C?té J, Laprise R. 2002. Spectral decomposition of two- dimensional atmospheric fields on limited-area domains using the discrete cosine transform (DCT) [J]. Mon. Wea. Rev., 130 (7): 1812–1829, doi:10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2.

        Evensen G. 1994. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics [J]. J. Geophys. Res., 99 (C5): 10143–10162, doi:10.1029/ 94JC00572.

        Evensen G. 2003. The ensemble Kalman filter: Theoretical formulation and practical implementation [J]. Ocean Dynamics, 53 (4): 343–367, doi:10.1007/s10236-003-0036-9.

        Evensen G, van Leeuwen P J. 2000. An ensemble Kalman smoother for nonlinear dynamics [J]. Mon. Wea. Rev., 128 (6): 1852–1867, doi:10.1175/1520-0493(2000)128<1852:AEKSFN>2.0.CO;2.

        Gaspari G, Cohn S E. 1999. Construction of correlation functions in two and three dimensions [J]. Quart. J. Roy. Meteor. Soc., 125 (554): 723–757, doi:10.1002/qj.49712555417.

        Hoke J E, Anthes R A. 1976. The initialization of numerical models by a dynamic-initialization technique [J]. Mon. Wea. Rev., 104 (12): 1551– 1556, doi:10.1175/1520-0493(1976)104<1551:TIONMB>2.0.CO;2.

        Kepert J D. 2009. Covariance localisation and balance in an ensemble Kalman filter [J]. Quart. J. Roy. Meteor. Soc., 135 (642): 1157–1176, doi:10.1002/qj.443.

        Lakshmivarahan S, Lewis J M. 2013. Nudging methods: A critical overview [M]// Park S K, Xu L. Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications (Vol. II). Springer: Berlin Heidelberg, 27–57, doi:10.1007/978-3-642-35088-7_2.

        Lei L L, Stauffer D R, Haupt S E, et al. 2012a. A hybrid nudging-ensemble Kalman filter approach to data assimilation. Part I: Application in the Lorenz system [J]. Tellus A, 64: 18484, doi:10.3402/tellusa.v64i0.18484.

        Lei L L, Stauffer D R, Deng A J. 2012b. A hybrid nudging-ensemble Kalman filter approach to data assimilation. Part II: Application in a shallow-water model [J]. Tellus A, 64: 18485, doi:10.3402/ tellusa.v64i0.18485.

        Lewis J M, Lakshmivarahan S, Dhall S. 2006. Dynamic Data Assimilation: A Least Squares Approach [M]. Cambridge, UK and New York, NY, USA: Cambridge University Press.

        Liu Y B, Bourgeois A, Warner T, et al. 2005. Implementation of observation-nudging based FDDA into WRF for supporting ATEC test operations [R]. 2005 WRF User Workshop. Paper 10.7.

        Liu Y B, Bourgeois A, Warner T, et al. 2006. An update on “obs-nudge”- based FDDA for WRF-ARW: Verification using OSSE and performance of real-time forecasts [R]. 2006 WRF User Workshop. Paper 4.7.

        Livings D. 2005. Aspects of the ensemble Kalman filter [D]. M. S. thesis, Reading University.

        閔錦忠, 王世璋, 陳杰, 等. 2012. 迭代EnSRF方案設(shè)計(jì)及在Lorenz96模式下的檢驗(yàn) [J]. 大氣科學(xué), 36 (5): 889–900. Min Jinzhong, Wang Shizhang, Chen Jie, et al. 2012. The implementation and test of iterative EnSRF with Lorenz96 model [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 36 (5): 889–900, doi:10.3878/j.issn.1006- 9895.2012.11185.

        Nair R D, Thomas S J, Loft R D. 2005. A discontinuous Galerkin global shallow water model [J]. Mon. Wea. Rev., 133 (4): 876–888, doi:10.1175/MWR2903.1.

        Ourmières Y, Brankart J M, Berline L, et al. 2006. Incremental analysis update implementation into a sequential ocean data assimilation system [J]. J. Atmos. Oceanic Technol., 23 (12): 1729–1744, doi:10.1175/ JTECH1947.1.

        Ruggiero G A, Ourmières Y, Cosme E, et al. 2015. Data assimilation experiments using the diffusive back-and-forth nudging for the NEMO ocean model [J]. Nonlin. Processes Geophys., 22 (2): 233–248, doi:10.5194/npg-22-233-2015.

        邵長(zhǎng)亮, 閔錦忠. 2015. 集合均方根濾波同化地面自動(dòng)站資料的技術(shù)研究 [J]. 大氣科學(xué), 39 (1): 1–11. Shao Changliang, Min Jinzhong. 2015. A study of the assimilation of surface automatic weather station data using the ensemble square root filter [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 39 (1): 1–11, doi:10.3878/j.issn.1006-9895. 1406.13263.

        Stauffer D R, Seaman N L. 1990. Use of four-dimensional data assimilation in a limited-area mesoscale model. Part I: Experiments with synoptic-scale data [J]. Mon. Wea. Rev., 118 (6): 1250–1277, doi:10.1175/1520-0493(1990)118<1250:UOFDDA>2.0.CO;2.

        Stauffer D R, Bao J W. 1993. Optimal determination of nudging coefficients using the adjoint equations [J]. Tellus A, 45 (5): 358–369, doi:10.1034/j.1600-0870.1993.t01-4-00003.x.

        王世璋, 閔錦忠, 陳杰, 等. 2013. 迭代集合平方根濾波在風(fēng)暴尺度資料同化中的應(yīng)用 [J]. 大氣科學(xué), 37 (3): 563–578. Wang Shizhang, Min Jinzhong, Chen Jie, et al. 2013. Application of iterative ensemble square-root filter in storm-scale data assimilation [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 37 (3): 563–578, doi:10.3878/j.issn. 1006-9895.2012.11186.

        Whitaker J S, Hamill T M. 2002. Ensemble data assimilation without perturbed observations [J]. Mon. Wea. Rev., 130 (7): 1913–1924, doi:10.1175/1520-0493(2002)130<1913:EDAWPO>2.0.CO;2.

        Yang S C, Corazza M, Carrassi A, et al. 2009. Comparison of local ensemble transform Kalman filter, 3DVAR, and 4DVAR in a quasigeostrophic model [J]. Mon. Wea. Rev., 137 (2): 693–709, doi:10.1175/2008MWR2396.1.

        Zhang F, Snyder C, Sun J Z. 2004. Impacts of initial estimate and observation availability on convective-scale data assimilation with an ensemble Kalman filter [J]. Mon. Wea. Rev., 132 (5): 1238–1253, doi:10.1175/1520-0493(2004)132<1238:IOIEAO>2.0.CO;2.

        Zou X, Navon I M, Ledimet F X. 1992. An optimal nudging data assimilation scheme using parameter estimation [J]. Quart. J. Roy. Meteor. Soc., 118 (508): 1163–1186, doi:10.1002/qj.49711850808.

        Implementation and Testing of a Hybrid Back and Forth Nudging Ensemble Kalman Filter (HBFNEnKF) Data Assimilation Method

        ZHU Haonan1, MIN Jinzhong1, 2, and DU Ningzhu1, 2

        1,,210044;2,,210044

        Based on the “Back and Forth Nudging” (BFN) and Ensemble Kalman Filter (EnKF) methods, a Hybrid BFN EnKF (HBFNEnKF) data assimilation method was designed and tested using a channel shallow water model and a global shallow water model, separately. Furthermore, the performances of the HBFNEnKF, Hybrid Nudging EnKF (HNEnKF), and Ensemble Square-Root Filter (EnSRF) methods are discussed, with model error considered. The results showed that the HBFNEnKF method retains the continuity and smoothness of HNEnKF, avoids the discontinuity and unbalance problem of EnSRF, and has the highest convergence speed. Through a single variable observation experiment, the advantage of HBFNEnKF was clear; that is, HBFNEnKF can maintain the balance between different model variables. A scale investigation on the increment field showed that, compared with EnSRF, HBFNEnKF produces a better assimilation result at larger scales, and avoids a number of spurious increments at medium and smaller scales.

        Data assimilation, Hybrid Back and Forth Nudging EnKF (HBFNEnKF), Ensemble Square-Root Fiter (EnSRF),Shallow water model

        1006-9895(2016)05-0995-14

        P413

        A

        10.3878/j.issn.1006-9895.1510.15214

        2015-06-17;網(wǎng)絡(luò)預(yù)出版日期2015-10-20

        朱浩楠,男,1991年出生,碩士研究生,主要從事中小尺度資料同化技術(shù)研究。E-mail: viczhn@gmail.com

        閔錦忠,E-mail: minjz@nuist.edu.cn

        國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃項(xiàng)目2013CB430102,國家自然科學(xué)基金重點(diǎn)項(xiàng)目41430427,國家自然科學(xué)基金項(xiàng)目41505089

        Funded by National Basic Research Program of China(Grant 2013CB430102), Key Program of National Natural Science Foundation of China (Grant 41430427), National Natural Science Foundation of China (Grant 41505089)

        猜你喜歡
        效果方法
        按摩效果確有理論依據(jù)
        學(xué)習(xí)方法
        迅速制造慢門虛化效果
        抓住“瞬間性”效果
        中華詩詞(2018年11期)2018-03-26 06:41:34
        可能是方法不對(duì)
        模擬百種唇妝效果
        Coco薇(2016年8期)2016-10-09 02:11:50
        用對(duì)方法才能瘦
        Coco薇(2016年2期)2016-03-22 02:42:52
        四大方法 教你不再“坐以待病”!
        Coco薇(2015年1期)2015-08-13 02:47:34
        賺錢方法
        捕魚
        蜜臀av 国内精品久久久| 国产午夜精品视频在线观看| 日韩女优精品一区二区三区| 亚洲精品国产suv一区88| 免费av片在线观看网站| 国产第一页屁屁影院| 就去吻亚洲精品欧美日韩在线| 在线播放国产女同闺蜜| 国产精品一区二区久久精品蜜臀| 国产免费人成视频网站在线18| 国产永久免费高清在线| 久久AV老司机精品网站导航| 国产成人精品亚洲午夜| 亚洲第一区二区快射影院| 中文字幕亚洲一区视频| 中文字幕一区二区中出后入 | 久久久久88色偷偷| 人人添人人澡人人澡人人人人| 毛片免费在线播放| 久久av一区二区三区下| 自拍视频在线观看首页国产| 久久久久久亚洲av无码蜜芽| 欧美激情二区| 女同中文字幕在线观看| 揄拍成人国产精品视频| 国产精品久久久久久无码| 精品久久久久久电影院| 特级黄色大片性久久久| 少妇人妻综合久久中文字幕| 伊人色综合视频一区二区三区 | 亚洲精品一区二区三区麻豆| 大地资源在线影视播放| 成人免费ā片在线观看| 扒下语文老师的丝袜美腿| 亚洲精品有码日本久久久| 中国丰满熟妇xxxx性| 国产伦精品一区二区三区四区| 国产精品专区一区二区av免费看| 国产视频一区二区三区在线免费| 欧美粗大猛烈老熟妇| 亚洲精品亚洲人成在线播放 |